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AN  LCAO  LOCAL  DENSITY  FUNCTIONAL  APPROACH  TO 
SURFACE  ELECTRONIC  STRUCTURE  CALCULATIONS 

By 

JOHN  WALLACE  MINTMIRE 
DECEMBER  I98O 
Chairman:   John  R.  Sabin 
Major  Department:   Physics 

A  new  method  for  calculating  the  electronic  structure  of  thin  films 
is  presented  using  the  local  density  functional  formalism  in  the  linear 
combination  of  atomic  orbitals  (LCAO)  approximation.   The  self-consistent 
solution  of  the  resulting  secular  equations  is  made  tractable  through  the 
use  of  fitting  procedures  for  approximating  the  charge  density  (for  the 
Coulomb  potential)  and  the  cube  root  of  the  charge  density  (for  the  ex- 
change potential)  as  linear  combinations  of  fitting  functions,  allowing 
the  construction  of  a  more  efficient  (in  terms  of  computational  effort) 
procedure  than  would  otherwise  be  possible.   Hermi te-Gauss ian  functions 
comprise  the  basis  sets,  with  different  sets  of  Hermi te-Gauss ian  functions 
used  for  the  orbital  basis  set,  for  the  charge  density  fitting  functions, 
and  for  the  functions  used  to  fit  the  cube  root  of  the  charge  density. 
Results  are  presented  for  two  systems,  the  atomic  hydrogen  monolayer  and 
the  atomic  beryllium  monolayer,  which  demonstrate  the  feasibility  and 


VI 


applicability  of  this  approach.   Techniques  are  discussed  for  improve- 
ment of  the  computational  procedure  in  order  to  enable  calculations  to 
be  made  on  larger  and  more  complex  systems  than  those  presented  in  this 
dissertation. 
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CHAPTER  I 
INTRODUCTION 


1-1.  General  Remarks 


The  study  of  the  electronic  properties  of  solid  surfaces  is  becoming 
Increasingly  more  important  as  techniques  dependent  on  surface  properties 
(such  as  catalysis  techniques  and  ultra-thin-film  devices)  occupy  more 
extensive  roles  in  modern  technology.   Experimental  interest  in  these 
systems  has  expanded  greatly  in  recent  years,  in  large  measure  due  to 
improved  vacuum  techniques  over  the  past  two  decades.   The  theoretical 
study  of  the  electronic  properties  of  surfaces  thus  faces  a  challenge  to 
develop  new  methods  for  understanding  the  behavior  of  these  systems. 

Many  experimental  methods  for  studying  surfaces  have  been  developed 
or  improved  in  the  last  decade  using  a  multitude  of  probe  techniques. 
Electron  beams  impinging  upon  the  surface  of  interest  are  used  in  Auger 
Electron  Spectroscopy  (AES) ,  Low  Energy  Electron  Diffraction  (LEED) ,  and 
the  Scanning  Electron  Microscope  (SEM)  to  yield  information  about  the 
electronic  structure  and  nuclear  geometry  of  the  surface,  v/hile  Electron 
Spectroscopy  for  Chemical  Analysis  (ESCA  or  XPS)  and  Ultraviolet  Photon 
Spectroscopy  (UPS)  use  photons  (X-ray  and  UV  respectively)  as  a  probe  to 
provide  additional  information  about  the  electronic  structure  of  surfaces 
(for  a  historical  review  of  these  methods,  see  Lee,  1977). 

in  this  work  we  propose  a  new  computational  approach  to  calculate 
the  approximate  electronic  two-dimensional  band  structure  and  total  energy 


of  a  slab  having  finite  thickness  in  one  direction  and  two-dimensional 
periodicity  in  the  plane  normal  to  the  direction  of  finite  thickness. 
The  slab  system  will  thus  possess  both  a  surface  and  the  periodicity  in 
directions  parallel  to  the  surface  characteristic  of  the  bulk  solid. 
Thus  by  increasing  the  thickness  of  the  slab,  the  properties  of  the  semi- 
infinite  solid  may  be  approached  as  a  limit. 

Methods  will  be  described  in  this  work  for  solving  the  secular 
equation  resulting  from  a  Linear  Combination  of  Atomic  Orbitals  (LCAO) 
tight-binding  approach  to  the  Xcx  formalism  (for  a  review  see  Connolly, 
1976).   The  computational  difficulties  arising  from  the  evaluation  of 
matrix  elements  involving  p^/^(r_)  (where  p  is  the  electronic  charge  den- 
sity) are    treated  by  fitting  p^'^  to  a  linear  combination  of  two- 
dimensional  ly  periodic  fitting  functions  in  a  manner  equivalent  to  that 
of  Dunlap  (Dunlap,  Connolly,  and  Sabin,  1979a  and  1979b)  for  molecular 
systems.   in  order  to  reduce  the  number  of  computationally  expensive 
integrals  which  must  be  evaluated,  the  electronic  charge  density  p  is  also 
fit  to  a  linear  combination  of  a  different  set  of  two-dimensional ly  per- 
iodic fitting  functions.   A  mathematical  formalism  is  derived  for  fitting 
p  based  on  the  procedure  introduced  by  Dunlap  (Dunlap,  Connolly,  and 
Sabin,  1979a)  and  Mintmire  (1979)  for  molecular  systems  which  finds  the 
variational  extremum  of  an  approximate  electron  repulsion  energy  func- 
tional of  the  fit  to  p.   This  approximate  electron  repulsion  energy 
functional  is  bounded  by  the  electron  repulsion  energy  of  p  with  itself 
and  is  equal  to  the  electron  repulsion  energy  of  P  vyith  itself  to  second 
order  in  the  difference  between  the  exact  density  P  and  the  fitted  density. 
The  formalism  derived  by  Dunlap  and  Mintmire  in  the  above  mentioned 
references  for  the  molecular  case  unfortunately  cannot  be  used  directly 


for  extended  systems  due  to  the  long  range  nature  of  the  Coulomb  inter- 
action resulting  in  certain  integrals  tending  toward  infinity  in  the 
limit  as  a  finite  molecular  system  increases  in  size  to  approach  the  two- 
dimensional  ly  periodic  slab.   Thus  changes  were  introduced  into  the 
original  formalism  to  yield  a  scheme  equivalent  to  the  original  molecular 
density  scheme  (Dunlap,  Connolly,  and  Sabin,  1979a),  but  with  proper 
1 imi  t  ing  behavior. 

This  new  computational  scheme  is  tested  on  two  systems:   the  atomic 
hydrogen  monolayer  and  the  atomic  beryllium  monolayer.   The  results  for 
the  hydrogen  monolayer  demonstrate  the  feasibility  of  the  method  and 
illustrate  some  of  the  effects  of  changes  in  basis  sets  and  in  the  number 
of  points  considered  in  the  surface  Brillouin  zone.   The  atomic  beryllium 
results  further  demonstrate  the  feasibility  of  the  method  and  selected 
results  are  compared  with  experimentally  determined  results  such  as  the 
work  function  and  Auger  spectra  of  the  beryllium  surface. 

1-2.  Historical  Background 

Early  calculations  of  the  electronic  structure  of  surfaces  (excluding 
cluster  calculation  approaches)  normally  fell  into  one  or  two  categories: 

1)  treating  the  system  as  an  electron  gas  and  neglecting  the 
electron-nuclear  lattice  interactions  (by  assuming  a  homogeneous 
positive  charge  density),  such  as  the  results  by  Lang  (1973) 
using  a  jellium  model  u'\th   a  step  potential  to  simulate  the 
potential  at  the  surface;  or 

2)  band  structure  calculations  taking  into  account  the  interactions 
between  electrons  and  the  nuclear  lattice,  although  these 


calculations  almost  universally  were  not  self-consistent  and 
frequently  were  semi-empi r ical  methods  using  experimental  data 
or  bulk  band  structure  data. 
We  wish  to  review  briefly  the  history  of  calculations  using  the  band 
structure  approach  (with  self-consistent  procedures)  over  the  last  several 
years  in  order  to  see  how  our  proposed  method  compares  with  previous 
work.   One  of  the  earliest  calculations  to  report  self-consistent  calcu- 
lations for  the  electronic  structure  of  a  thin  film  system  was  by 
Alldredge  and  Kleinman  (1972).   These  two  workers  calculated  the  band 
structure  of  a  13  atom  thick  layer  lithium  system,  taking  into  account 
the  effects  of  the  nuclear  lattice,  using  a  non-self-consistent  plane 
wave  solution  which  was  later  extended  to  a  self-consistent  scheme  for 
calculations  on  the  same  lithium  system  (Alldredge  and  Kleinman,  1974). 
The  lack  of  translat ional  symmetry  in  the  direction  normal  to  the  surface 
was  treated  by  forcing  all  orbitals  to  vanish  at  planes  located  some 
large  distance  from  the  surface  (Alldredge  and  Kleinman  used  planes  lo- 
cated three  times  the  layer  separation  away  from  the  exterior  nuclear 
site.),  equivalent  to  invoking  an  impenetrable  barrier  for  electrons 
away  from  the  surface.   Core  states  were  not  treated  directly  in  this 
method,  but  core  effects  were  included  in  a  pseudopotent ial  model. 

Applebaum  and  Hamann  (1972)  introduced  a  model  of  the  semi  - inf i ni te 
surface  by  requiring  all  orbitals  to  match  the  bulk  orbitals  as  boundary 
conditions  on  one  side  of  a  thin  film  system.   Model  potentials  including 
core  effects  were  used  to  generate  orbitals  in  the  Laue  representation 
(von  Laue,  1931),  where  a  Fourier  expansion  in  the  parallel  coordinates 
is  combined  with  the  coordinate  representation  in  the  normal  direction, 
using  a  numerical  grid  for  the  coordinate  representation  component  of 


all  evaluated  functions.   This  approach  was  used  primarily  to  study  re- 
laxation effects  on  the  (ill)  surface  (Appelbaum  and  Hamann,  1973)  and 
the  (100)  surface  (Appelbaum,  Baraff,  and  Hamann,  1975a  and  1975b)  of 
s  i 1 i  con . 

Cooper  (1973)  reported  calculations  on  the  copper  (100)  monolayer 
using  a  muffin-tin  version  of  the  KKR  method  adapted  to  two-dimensional ly 
periodic  systems.   The  lack  of  translational  symmetry  in  the  surface 
normal  direction  was  treated  with  an  impenetrable  barrier  in  a  similar 
manner  to  that  of  Al Idredge  and  Kleinman  (1972).   Various  refinements  to 
this  basic  approach  were  introduced  by  Kohn  (1975)  and  Kar  and  Soven 
(1975)  who  simultaneously  reported  muffin-tin  KKR  approaches  that  did  not 
require  the  use  of  an  impenetrable  barrier  located  away  from  the  surface. 
Krakauer  and  Cooper  (1977)  further  adapted  this  approach  into  a  local 
combination  of  muffin-tin  orbitals  (LCMTO)  approach  incorporating  non- 
muffin-tin  corrections  to  the  potential. 

Gay,  Smith,  and  Arlinghaus  (1977)  introduced  a  sel f -cons i stent  band 
structure  method  in  an  LCAO  approach  using  Gaussian  basis  functions. 
Their  approach  utilized  an  effective  potential  constructed  from  two  terms: 

1)  the  initial  potential  generated  by  overlapping  atomic  potentials; 
and 

2)  differences  between  the  effective  potential  at  each  iteration 
in  the  self-consistent  process  and  the  initial  potential  cal- 
culated as  a  Fourier  expansion. 

This  approach  has  been  very  successful  in  describing  the  band  structure 
of  transition  metals,  leading  to  computed  results  for  a  nine  layer  model 
of  the  copper  (100)  surface  (Smith,  Gay,  and  Arlinghaus,  I98O)  and  a 
nine  layer  model  of  the  nickel  (100)  surface  (Arlinghaus,  Gay  and  Smith, 


1980).   This  method  unfortunately  appears  not  to  admit  a  straightforward 
method  for  evaluating  the  total  energy  of  the  system  due  to  the  nature 
of  the  Fourier  transform  techniques. 

Wang  and  Freeman  (1978)  introduced  another  LCAO  approach  for  thin 
films  based  on  the  Discrete  Variation  Method  (DVM)  of  Ellis  and  Painter 
(1971)  where  the  orbitals  are  expressed  in  a  numerical  basis  set.   The 
charge  density  used  for  generating  the  secular  matrix  of  each  iteration 
in  this  self-consistent  procedure  is  constructed  by  fitting  the  current 
true  charge  density  to  a  charge  density  resulting  from  superimposed  over- 
lapping spherically  symmetric  atomic  charge  densities.   This  method  has 
been  used  to  calculate  the  band  structure  of  one,  three,  and  five  layer 
models  of  the  nickel  (001)  surface  (Wang  and  Freeman,  1979)  and  a  nine 
layer  model  of  the  nickel  (001)  surface  (Wang  and  Freeman,  I98O).   No 
total  energies  have  been  reported  using  this  method,  although  no  inherent 
difficulties  appear  to  prevent  this  approach  from  yielding  total  energies 
in  a  straightforward  fashion. 

In  addition  to  the  method  introduced  by  them  mentioned  previously, 
an  alternative  self-consistent  LCAO  procedure  using  Gaussian  basis  func- 
tions has  been  introduced  by  Appelbaum  and  Hamann  (1978),  with  calculated 
results  for  the  copper  (100)  monolayer.   Their  approach  is  similar  to  that 
derived  in  this  work  in  the  use  of  fitting  techniques  for  the  charge  den- 
sity p  and  for  p^'^.   Appelbaum  and  Hamann  fit  the  charge  density  with  a 
least-squares  fitting  procedure  to  a  periodic  sum  of  spherically  symmetric 
Gaussians  using  "floating"  Gaussians  at  sites  other  than  the  nuclear 
centers  in  order  to  describe  asymetric  components  of  the  charge  density. 
Using  the  same  basis  of  spherically  symmetric  Gaussians  as  used  for 
fitting  the  charge  density,  Appelbaum  and  Hamann  then  fit  the  sum  of  the 


local  density  functional  exchange  potential  and  the  total  Coulomb  poten- 
tial minus  a  screened  nuclear  attraction  term,  which  is  treated  outside 
the  fitting  procedure.   This  approach  has  been  used  to  calculate  the 
band  structure  of  the  eleven  layer  model  of  the  titanium  (0001)  surface 
(Feibelman,  Appelbaum,  and  Hamann,  1979)  and  the  same  titanium  system 
with  a  chemisorbed  hydrogen  monolayer  (Feibelman  and  Hamann,  I98O)  and 
with  a  chemisorbed  nitrogen  monolayer  (Feibelman  and  Himpsel,  I98O). 
Recent  calculations  (Feibelman,  Hamann,  and  Himpsel,  I98O)  report  cohesive 
energy  curves  for  the  eleven  layer  model  of  the  titanium  (OOOl )  surface 
with  a  chemisorbed  hydrogen  monolayer  at  various  distances  from  the  ex- 
terior layer  of  titanium.   Unfortunately  complete  details  as  to  how  the 
total  energy  terms  were  evaluated  in  this  formalism  have  not  been  pub- 
I i  shed  at  this  time. 

The  formalism  described  in  this  work  is  similar  to  the  techniques 
discussed  by  Feibelman,  Appelbaum,  and  Hamann  (1979);  the  primary  dif- 
ferences arise  from  the  fitting  techniques  utilized.   In  our  approach 
a  basis  set  of  Gaussian  functions  (not  restricted  to  spherically  symmetric 
Gaussians)  is  used  to  fit  the  p^'  ^  portion  of  the  local  density  functional 
exchange  potential  by  using  least-squares  fitting  methods.   The  charge 
density  is  then  fitted  using  the  variational  methods  mentioned  in  Section 
1-1,  which  results  in  the  evaluated  Coulombic  energy  portion  of  the  total 
energy  agreeing  with  the  Coulombic  energy  resulting  from  the  non-fitted 
density  to  second  order  in  the  error  in  the  charge  density  introduced  by 
the  fitting  procedure.   We  feel  this  is  a  more  rational  scheme  than  a 
least-squares  fitting  method  for  fitting  the  charge  density  when  the 
fitted  charge  density  will  be  used  to  construct  the  Coulombic  potential 
due  to  the  electron  charge  and  to  evaluate  the  Coulombic  contribution 
to  the  total  energy. 


1-3.  Organization  of  Dissertation 

There  are  two  fundamental  topics  to  cover  in  this  dissertation. 
The  first  topic  is  a  basic  description  of  the  electronic  band  structure 
approach.   The  mathematical  formalism  of  our  approach  is  presented  in 
Chapter  11,  followed  by  a  brief  description  of  the  computational  tech- 
niques used  to  implement  the  formalism  given  in  Chapter  III.   The  results 
are  the  second  fundamental  topic  to  be  discussed.   Computational  results 
for  the  atomic  hydrogen  monolayer  and  the  atomic  beryllium  are  presented 
and  discussed  in  Chapter  IV,  with  a  summary  of  our  method  and  results 
in  Chapter  V. 

The  appendices  at  the  end  of  this  dissertation  discuss  in  extended 
detail  several  topics  relating  to  the  computational  procedure.   Appendix 
A  defines  the  Hermi te-Gauss i an  functions  used  in  the  work  and  integrals 
involving  these  functions  are  reduced  to  forms  which  may  be  evaluated 
using  algorithms  described  in  this  appendix  and  later  appendices.   The 
multipole  expansion  techniques  used  In  evaluating  Coulomb  Integrals  are 
discussed  In  Appendix  B,  while  equations  necessary  for  the  evaluation  of 
integrals  using  reciprocal  lattice  expansions  are  derived  in  Appendix  C. 
Finally,  Appendix  D  describes  the  three-dimensional  grid  used  for  the 
numerical  integration  algorithms  which  are  part  of  the  computational 
fitting  procedure. 


CHAPTER  I  I 
MATHEMATICAL  FORMALISM  AND  THEORY 


2-1.  Symmetry  Properties  of  Thin  Films 

As  mentioned  in  Chapter  I,  the  physical  systems  of  actual  interest 
were  studied  using  an  ideal  periodic  thin  film  model.   Given  an  ideal 
three-dimensional  periodic  lattice  of  nuclei,  the  thin  film  may  be  vis- 
ualized as  that  portion  of  the  three-dimensional  lattice  enclosed  between 
two  parallel  planes  of  finite  separation  as  illustrated  in  Figure  2-1. 
A  more  formal  definition  of  the  systems  under  consideration  is  expressed 
by  the  following  two  conditions  on  the  positions  of  the  nuclei  In  the 
nuclear  lattice: 

1)  There  exist  two  Independent  vectors  R^.  and  R^„  such  that  a 
translation  of  coordinates  by  R_  =  mR^.  +  n R^„ ,  where  m  and  n  are 
integers,  results  In  an  equivalent  lattice; 

2)  and  there  exist  two  planar  surfaces,  parallel  to  the  plane  of 
periodicity  containing  R,.  and  R^^ ,  with  finite  separation  such 
that  all  of  the  nuclear  sites  are  contained  within  the  region 
enclosed  by  the  two  parallel  planes. 

The  first  condition  requires  two-dimensional  periodicity  while  the  second 
condition  requires  that  the  thin  film  indeed  has  finite  thickness. 

The  set  of  all  possible  translations  R^, 

R^  =  mR^.  +  nR_„  ,   (m  and  n  integers)  (2-1) 
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Figure  2-1.  Schematic  representation  of  a  side  view  of  a  two-d i mens  ion- 
ally  periodic  slab  v/ith  three  layers. 
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defines  a  two-dimensional  geometric  lattice  of  points,  the  planar  Bravais 
lattice.   There  are  five  possible  types  of  planar  Bravais  lattices,  which 
have  different  point  group  symmetries  associated  with  them  in  addition  to 
the  translat ional  symmetry,  as  illustrated  in  Figure  2-2.   An  elementary 
treatment  of  these  systems  is  given  in  Blakemore  (197^)- 

Once  a  planar  Bravais  lattice  is  chosen,  the  treatment  of  symmetry 
is  similar  to  that  of  the  three-dimensional  crystalline  lattice  as  des- 
cribed by  Tinkham  (1964)  or  by  Slater  (1972a).  A  unit  cell  may  be  con- 
structed as  shown  in  Figure  2-3  in  a  manner  analogous  to  that  of  the 
three-dimensional  unit  cell.   Since  the  thin  film  is  actually  a  three- 
dimensional  system  while  the  planar  Bravais  lattice  is  only  a  two- 
dimensional  geometrical  construct,  the  thin  film  unit  cell  is  a  parallel- 
opiped  which  is  finite  in  directions  parallel  to  the  plane  of  periodicity 
but  infinite  in  extent  in  the  directions  normal  to  the  plane  of  period- 
icity.  The  arrangement  of  nuclear  centers  inside  each  unit  cell  is 
unrestricted  except  for  the  previous  constraint  that  the  thin  film  have 
finite  thickness.   Therefore  the  complete  arrangement  of  nuclei  need  not 
have  the  full  point  group  symmetry  of  the  Bravais  lattice. 

Thus  the  symmetry  group  of  the  total  system  (or  rather  the  plane 
group,  since  we  are  currently  considering  the  set  of  two-dimensional 
symmorphic  symmetry  operations)  is  a  subgroup  of  the  plane  group  of  the 
Bravais  lattice.   Using  an  approach  similar  to  that  of  Koster  (1957)  the 
plane  group  consists  of  symmetry  operators  of  the  form  {P|R^}  which  cor- 
respond to  a  point  group  coordinate  transformation  followed  by  a  trans- 
lation of  coordinates: 

{pIr}  r  =  Pr  +  R  (2-2) 
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igure   2-2.    The    five    tv/o-d imens ional    lattices   and    their   symmetries. 

(From  Sol i  d  State  Phys  ics  by  J.  S.  Blakemore.  Copyright  (c) 
1969  by  U.  B„  Saunders  Company.  Reprinted  by  permission  of 
Holt,    Rinehart,    and   Winstono) 
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Figure  2-3.  1  I  1 ust rat  ion  of  parall elopiped  unit  cell  for  the  two-dimen- 
sional ly  periodic  monolayero   Note  that  the  unit  cell  extends 
to  infinite  distance  in  both  directions  parallel  to  the  z 
axi  s. 


The  operator  {P|R^}  can  be  seen  to  have  the  inverse  operator  {p|r} 
such  that: 

{pIr}"'  =  {P"'|-P"^  R.)  (2-3) 

Given  these  definitions  the  operation  of  a  symmetry  operator  upon  a 
function  f ( 0  may  be  defined: 

{P|R}  f(r_)  =  f({P|R}"'  r)  {2-k) 

Bloch's  theorem  applies  to  the  two-d imens iona lly  periodic  thin  film 
in  a  manner  analogous  to  that  for  the  three-dimensional  crystalline  case. 
From  group  theory  (Tinkham,  196'*)  one  can  show  that  for  an  operator  M 
which  commutes  with  all  the  translation  operators  R  (where  R  denotes  the 
subgroup  of  translation  operators  of  the  plane  group  given  by  {e|r^},  if 
E  denotes  the  identity  point  group  operator)  the  eigenfunct ions  of  M  can 
be  chosen  as  members  of  the  irreducible  representations  of  the  translation 
group.   These  irreducible  representations  may  be  labelled  with  wave 
vectors  j<,,  that  lie  parallel  to  the  plane  of  periodicity  such  that  if 
the  function  <l>(r_;  Ji//)  belongs  to  the  J<-,  irreducible  representation, 
then: 


R  ^in    k^p    =   ^{r_  -    R;  k^^)  =  exp(ik//-R)  $(_r;  k^/)  (2-5) 

In  a  manner  similar  to  the  three-dimensional  crystalline  case,  the 
set  of  wave  vectors  k_.  .    is  not  a  set  of  un  ique  labels  for  each  individual 
irreducible  representation.   For  example,  if  the  vector  c_   is  normal  to 
both  primitive  lattice  vectors  R^.  andR^  (e.g.,£=R^.  xR^),  then  two 
reciprocal  lattice  vectors  may  be  defined  as: 
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-»  "  R 1  •  ( R2  X  c )  ^^2  ^  ^^  (2-6) 


and 


-2  "  R2-(R^  X  c)   ^^1  '^  £)  (2-7) 

These  reciprocal  lattice  vectors  have  the  property  that: 


R.-K.  =  27r  6.  . 
-'  -J       ij 


(2-8) 


Thus  if  we  define  a  general  reciprocal  lattice  vector  K   such  that: 

K  =  hKj  +  kK^  (2-9) 

then  any  two  wave  vectors  k..   and  k'       related  by: 

h.'//   =  ii//  +  Ji  (2-10) 

label  the  same  irreducible  representation  since: 

$(£  +  R;  k^p    =  exp[i(k^/  +  K).R]  $(^;  k^^)  (2_n) 

=  exp[il<^^.R^]  $(£;  k   ) 

Thus  we  may  restrict  our  choice  of  wave  vectors  to  the  two-dimensional 
central  Brillouin  zone  of  wave  vectors  k^^  closer  to  the  origin  than  to 
any  other  point  described  by  an  arbitrary  reciprocal  lattice  vector  K  ?^  0 . 

With  the  Brillouin  zone  defined  we  may  restate  Bloch's  theorem 
in  a  more  straightforward  form  for  use  in  solving  finite  secular  matrix 
problems.   The  matrix  elements  of  an  operator  M,  which  commutes  with  all 
translation  operators  R,  will  be  zero  for  matrix  elements  between  functions 
belonging  to  different  irreducible  representations  labelled  by  k   and  k'  . 
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Recalling  that  the  symmetry  operators  are  unitary,  or: 

<$.(£_)  I  r"'  ^Ar)>  =   <R  $j(£)|  $.(£)>  (2-12) 

We  next  consider  the  matrix  element 


^-^i^L''  ii//)l  M  i.W;   !<//)>  =  <$.(n  k^/)|  r"'  r  m  $.(n  k^^)>  (2-13) 

=  <R  (I>.(n  k^^)|  M  R  *.(£;  k^^)> 

=  exp[i(k^/  -  k//)-R]  <*j(l;  k^/)|  M  $(£;  k^^)> 

Since  R^  may  represent  any  lattice  point  and  by  the  hypothesis  we  are 
given  that  j<^.  -  j<,,  ?^  j<,  then  the  matrix  element  must  equal  zero: 

<*,(£;  Ji//)|  M  <1'.(l;  Ji//)>  =  0  (2-Ii») 

A  star  of  the  wave  vector  k_..   may  be  defined  as  the  set  of  wave 
vectors  generated  by  the  operation  of  each  symmetry  operator  in  the  point 
symmetry  subgroup  of  the  plane  group  of  the  thin  film  system  on  the  wave 
vector  Ji//-   Since  the  wave  vectors  k_,,  are  merely  labels  for  the  Irre- 
ducible representations  of  the  translation  group,  it  is  understood  that 
if  a  symmetry  operation  yields  a  wave  vector  Wl.   outside  the  central 
Brillouin  zone  then  the  v;ave  vector  k'..   =  Wl,   +   K  equivalent  to  the  wave 
vector  J<'1,  is  used  as  the  member  of  the  star.   Thus  it  can  be  shown  that 
the  matrix  elements  of  an  operator  H  which  commutes  with  the  point  group 
symmetry  operators  P  (where  P  denotes  {P|0_})  have  the  property  that: 

<*•,(£;  ii//)|  M  $.(l;  k//)>  =  «J'j(ji;  k^/)|  P"'  P  M  ^.{r_;    k^^)>  (2-15) 

=  <P$.(_r;  k^^)\    M  P  i> .  in    k//)> 


17 

V/e  notice  that  the  new  function  P'? .  (£;  V.^)    belongs  to  the  irreducible 
representation  labelled  by  PJ<//  ,  since  if  R  is  a  general  translation 
operator: 

R  P  $.(r;  k//)  =  R  $.(p-l  r;  k^^)  (2_,6) 

=  $.(P"'  [r  -  R];  k^^) 

=  exp[-ik^^-P"'  R]  t.(p"^  £;  k   ) 

=  exp[-i(Pk^^)-R]  P  .|>.(_r;  k^^) 

From  group  theory  we  recall  that  the  eigenvectors  of  an  operator  M 
which  commute  with  the  operators  of  the  symmetry  group  will  be  degenerate 
if  a  wave  vector  star  possesses  more  than  one  unique  member,  since  each 
irreducible  representation  V^^  in  the  wave  vector  star  will  contain  one 
degenerate  eigenf unct ion .  Equation  2-16  implies  that  this  relationship 
is  also  true  for  the  finite  secular  matrix  problem  if  the  basis  sets  are 


such  that  if  the  basis  set  for  irreducible  k,,  contains  a  functi 

Pi// 


j<,,  contains  a  function 
*;(i;  ii//).  then  the  basis  set  for  the  Pl^  .  irreducible  representation 


must  contain  the  function  P^. (£;  k.,). 


2-2.  The  Xa  Method 


The  Xa  method  is  an  approximate  method  for  calculating  the  electronic 
structure  of  atoms,  molecules,  and  extended  systems  using  an  approach 
originated  by  Slater  (1951).   The  Born-Oppenheimer  approximation  (192/) 
is  used  to  simplify  the  problem  of  finding  the  eigenvalues  of  the  elec- 
tronic Hamiltonian: 
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[H(R_)  -  E(R)]y(_r.)  =  0  (2-17) 

where  (using  Hartree  atomic  units  throughout): 

S(R)=z(-ivJ.  I  i  .  r|^|.  z^,        (2-,8) 

I     —I   j>i   ij    m'— 1  — n '   m>n'-m-fi' 

The  coordinates  R   refer  to  nuclear  coordinates  and  the  coordinates  r. 
— m  _i 

refer  to  electronic  coordinates. 

The  Xa  method  originated  as  a  means  of  approximating  the  Hartree- 
Fock  model  by  replacing  the  non-local  exchange  terms  with  a  local  func- 
tional of  the  electron  density  (the  origins  of  the  Xa  method  are   discussed 
at  length  in  Slater  (IS?'^)).   A  different  interpretation  of  the  local 
exchange  potential  may  be  made  using  an  approach  based  on  the  local  density 
functional  formalism  of  Kohn  and  Sham  (I965).   This  local  density  func- 
tional formalism  is  based  on  the  theorem  by  Hohenberg  and  Kohn  (196'»). 
The  Hohenberg-Kohn  theorem  states  that  for  a  Hamiltonian  of  the  general 
form 

H  =  Z{-i  vj  +  v(£.)  +      I     -^  }  (2-19) 

i     -i  j>i   ij 

there  exists  an  energy  functional  G[p]  such  that  if  p{r)    is  the  electron 
density  resulting  from  Y(r.),  then 


—I 


.3n 


G[p]  =  /d^"r.  y"(£.)  [H  -  vv(£.)]  y(j^.)  (2-20) 


I    —I       .  — I 
I 


Furthermore  if  Poir)    corresponds  to  the  ground  state  'fo(£-)  of    the  Ham- 
iltonian H,  then  the  total  energy  functional: 

E[p]  =  G[p]  +  /d^r  p(£)  v(£)  (2-21) 

has  a  variational  minimum  about  the  electron  density  po(r). 
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If  instead  of  considering  just  the  electron  density  0(0  but  considering 
the  first  order  density  matrix  P  (_r ;  r_' )  which  may  without  loss  of  gener- 
ality be  assumed  to  be  of  the  form: 

p(_r;  £')  =  Z  n.  $?(£')  $.(£)  (2-22) 

i 

where  the  functions  $.  form  an  orthonormal  basis  and  the  occupation  num- 

I 

bers  have  the  property  0  <.  n.  <.  1,  then  we  may  define  a  new  energy  func- 
tional E   [p]  : 
xc 

E^Jp]  =  G[p]  -  !d\   jd\'    6(r  -  r')[-jvj]  p  ( ^  r ' ) 

,   ,    ,    p(r)  p(r') 
4  /d^  jd^r<   = = (2-23) 


where  p  ( r_)  =  p  ( r_;  r)  .      The  second  term  on  the  right-hand  side  of  Equation 
2-23  corresponds  to  the  expectation  value  of  the  kinetic  energy  of  the 
first  order  density  matrix  p  (_r ;  r_')    and  the  third  term  corresponds  to  the 
Coulomb  interaction  of  the  charge  density  p (0  with  itself.   Thus  we  see 
that  the  energy  functional  E   [p]  will  correspond  to  the  exchange  and 
correlation  energy  components  of  the  total  energy. 

The  local  density  functional  formalism  is  based  on  the  approximation 
of  the  exchange  correlation  energy  functional,  which  has  an  unknown 
analytic  form,  by  a  prescribed  functional  of  the  form: 

E   [p]  -  Jd^r  {pt(r)  V   (pf)  +  pl(r)  V   (pO)  (2-2^*) 

xc      -^        —   xc         —   xc 

where  V   (r_)  is  a  function  of  the  electron  density  of  spin  up  for  pt  and 
of  the  electron  density  of  spin  down  for  p4-.   The  potential  V  (r)    is 
most  properly  denoted  as  the  exchange-correlation  potential.   However, 
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for  brevity  and  historical  continuity  we  shall  refer  to  V   (r)  henceforth 

xc  — 

as  the  exchange  potential. 

In  the  Xa  model  the  exchange  potential  V   is  chosen  to  be  of  the 
form  introduced  by  Slater  (1951)  with  the  addition  of  an  adjustable 
parameter  a   with  a  value  between  2/3  and  1  (Slater  and  Wood,  1971): 

1/3 


\M^=-i^^-^^  (2-25) 


For  convenience  we  shall  henceforth  refer  to  an  exchange  potential  operator 

V   which  when  operating  on  p(r)  has  the  form: 
xa  r       3     K  _ 

V  p  =  V^  (p+)p+  +  V   (p4-)p+  (2-26) 

xa     Xa         xa 

This  set  of  approximations  yields  an  energy  functional  of  the  form: 

}         7.          T.                                    9            pir)    p{r_') 
E[p]  =  J   /d\  /d^r'  {6(£  -  £')  V^  p(£;  r ' )  +  1   .  ^,| 


^^  /^'^  ^^xa  -  '  yr^Tj'  ^^r)^    ^  ^^      (2-27) 

m  '—   -ffl '         m>n  '  m    n  ' 

where  the  total  energy  includes  the  nuclear  repulsion  term  introduced  by 
the  Born-Oppenheimer  approximation. 

The  orbitals  $. (r)  are  chosen  by  assuming  as  an  approximate  form  the 
linear  combination  of  a  finite  set  of  basis  functions: 


<J>.(r)  =  Z  c  .  v.(r)  (2-28) 

I  —       mi   I  — 
m 

and  the  energy  functional  E[p]  is  minimized  using  the  variational  principle 
subject  to  the  constraint  that  the  orbitals  <J.  are  orthonormal  ized.   This 
procedure  yields  the  familiar  secular  equation: 

EH   c.=ES   c.e.  (2-29) 

mn   n I      mn   n i   i 
n  n 
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where  £.  are  the  orbital  eigenvalues  and 


S   =  /d^r  /(r)  .;  (r)  (2-30) 

mn   •'     m  —   n  — 

mn   ■'     m  —   err  —   n  — 

The  one-electron  Hamiltonian  H  rr('')  'S  given  by: 

err  — 

m  '—   — m '  '—   —  ' 

where  p  corresponds  to  p+  If  li*  and  v  correspond  to  spin  up,  and  p 
o  m      n  a 

corresponds  to  p4-  if  ii     and  0  correspond  to  spin  down. 

m      n 

2-3.  LCAO  Methods  for  Extended  Systems 

For  infinitely  extended  systems  the  limit  of  the  total  energy  of  a 
finite  section  of  the  infinite  system  as  the  size  of  the  finite  section 
is  systematically  increased  is  not  a  useful  definition  of  the  total  energy 
of  the  periodic  system.   It  is  easily  seen  that  the  total  energy  defined 
in  this  fashion  will  in  most  cases  increase  linearly  with  the  cross- 
sectional  area  of  the  finite  section  (for  two-d imens iona 1 ly  periodic 
systems  as  defined  in  Section  2-1)  and  thus  will  not  have  a  convergent 
limit.   The  total  energy  per  unit  cell  is  a  far  more  useful  quantity  that 
is  usually  convergent  in  the  above  mentioned  limiting  process  for  systems 
that  possess  zero  total  charge.   In  the  Xa  model  the  functional  yielding 
the  energy  per  unit  cell  may  be  expressed: 
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E[p]  =  -J   /^d^r  /d^r'  6(j^  -  £-)  vj  p  ( n  j^- ) 

1  ,      o       ,    o  P(r)  P(r')  Z 

+  1  Ld\  d\'    —^ =:r--  E/d^r   (r)  E  -, OL 

2  ■^  fi  r  -  r '     „  ''n    p  '-'     r  -  R   -  R 

—   —      R^  m  '—   -m   — ' 

*//^Vt'4^   ^'  |R   -"r"  -  R|  (2-33) 

R  m,n  '— m   -n   — ' 


where  the  primed  summation  over  the  indices  m  and  n  indicates  that  the 
terms  corresponding  to  m  =  n  when  R^  =  0  are  to  be  omitted.   The  inte- 
grations over  the  region  H  denote  an  integration  over  one  three-dimensional 
unit  cell  (as  defined  in  Section  2-1). 

For  computational  convenience  we  assume  Born-von  Karman   (1912) 
periodic  boundary  conditions,  such  that  all  eigenf unct ions  of  periodic 
operators  are  periodic  over  translations  of  the  general  form  M  R^  +  N  R^  , 
where  M  and  N  are  integers  and  R^.  and  R^„  are  the  primitive  lattice  vectors 
for  the  system  of  interest.   Periodic  boundary  conditions  reduce  the 
problem  involving  an  infinite  number  of  electrons  to  the  problem  of  a 
finite  number  of  electrons  contained  in  M  by  N  unit  cells,  such  that  the 
orbitals  obey  periodic  boundary  conditions.   This  yields  a  procedure 
which  is  effectively  a  cluster  calculation  on  an  M  by  N  unit  cell  system 

« 

with  an  external  environment  outside  being  a  periodically  repeated  version 
of  the  internal  system  rather  than  empty  space. 

Let  us  approximate  the  orbitals  '5.(r)  by  a  linear  combination  of  a 
finite  number  of  orbitals  t;^  (_r )  .   The  spin-up  one-electron  effective 
Hamiltonian  for  this  procedure  may  be  expressed, 

H^ff  (L)  =  T  ^'  4  /^'^'  F^=;-T  ■  I   '    |r-R  -R|  -  ^"t-if  ^ 

'—   —  '     _R  m  '—  — m  — ' 

(2-3^ 
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with  the  spin  down  effective  Hami 1  ton ian  having  an  analogous  form.   From 
equation  2-1^  we  see  that  the  logical  choice  of  basis  functions  will  be 
Bloch  functions  ^    (rj;   Ji//)  .  or  now  labelling  functions  with  the  wave 
vector  }<_,/■ 

p(_r;    r_')    =   l   I  n;(li//)    ^^in   Ji//)    ^i^L'   Ji//)  (2-35) 

i   ii// 

^(nii//)   =^   c,^(k//)    ijn    k//)  (2-36) 

m 

We  see  that  the  Born-von  Karman   periodic  boundary  conditions  require  that 

the  orthonormal izat ion  condition  be  expressed: 

where  the  integration  region  Q'    corresponds  to  the  region  defined  by  the 
M  by  N  unit  cells.   The  periodic  boundary  conditions  also  restrict  the 
choice  of  wave  vectors  to  a  discrete  set,  since 

.];.(£+  MRj  +  NR^;  k//)  =  exp[ik^^-(MRj  +  NR^)  ]  ^.W;    k^/)     (2-38) 

=  i^.in   k//) 

This  impl ies 

Ji//  "  M  Ji]  "^  '^-2'    ^"^   ^"*^  "  integers)  (2-29) 


are  the  only  allowed  choices  of  k_//,  in  addition  to  the  condition  that  we 
only  consider  wave  vectors  in  the  central  Brillouin  zone. 

The  Bloch  basis  functions  'ii.ir;   }<_,,)    may  be  constructed  from  localized 
f  unct  ions  U  .  (_r)  : 

^;(l;  ii//)  =  ^  exp[ik^^-R^]  U.(_r  -  R_)  (2-40) 

R 
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This  yields  a  secular  equation  of  the  form 

^  H^^(k//)  c^i(l<//)  =  .,  Z  S^^(k//)  c^.(l<//)  (2-Z.l) 

n  n 

where 

S^n^ii//^  =  ^  exp[ik^/-R]  /d^r  U_^(r+R)  U^  (r_)  (2-A2) 


H^in^i//)  =  ^  exp[ik//.R]  /d^r  U^  ( r_+R )  H^^^(0  U^(0  (2-43) 

R 


2-4.  Fitting  Methods  for  Thin  Films 


The  effective  Hamiltonian  in  Equation  2-4l  contains  a  term  propor- 

1/3 
tional  to  p   ;  one  cannot  usually  express  integrals  involving  this  term 

in  an  analytic  closed  form.   One  way  of  alleviating  this  difficulty  for 

1/3 
molecular  systems  is  the  approximation  of  the  functions  p    (r)  by  a  linear 

combination  of  fitting  functions  G  (r) ,  henceforth  referred  to  as  the 

m   — 

exchange  fit, 

p^/^r)  -Eg   G  (r)  {2-hk) 

—       m  m  — 
m 

This  approach  to  solving  the  problems  involving  the  analytic  behavior  of 

1/3 
p    was  originally  suggested  by  Sambe  and  Felton  (1975)  with  extensions 

to  the  formalism  and  computational  method  later  introduced  by  Dunlap, 
Connolly,  and  Sabin  (1977)  and  Mintmire  (1979).   In  these  previous  works 
the  total  charge  density  p  ( r_)  was  also  approximated  using  a  linear  com- 
bination of  a  set  of  charge  fitting  functions  F  (r),  this  procedure  being 
henceforth  referred  to  as  the  charge  fit, 

p(j:)  =^  ^  f„  ~^Jl)  (2-45) 
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The  rationale  behind  approximating  p  (_r )  is  to  reduce  the  number  of  analytic 
integrals  which  must  be  computed,  thus  achieving  a  savings  in  computa- 
tional effort. 

For  extended  systems  the  fitting  functions  G  (r)  and  F  (r)  are  com- 

m  —       m  — 

posed  of  periodic  sums  of  localized  functions  G  (r)  and  F  (r),  (where  a 

m  —       m  — 

localized  function  has  the  property  that  lim  r-^  G  (r)  =  0),  that  may  be 

r->TO     m  — 

expressed: 

G  (r)  =  V  G  (r  -  R)  {l-kG) 

m-    ^  m-   - 


m-    R  m-   - 


The  exchange  fitting  coefficients  g  are   chosen  using  a  least  squares 
fitting  procedure  in  which  the  set  of  conditions 


^  [  Ld^{p^/^(r)  -Eg   G  (r)}2]  =  0  (2-/*8) 

'  \/  ^—  m   m  — 


3g.   •'fi    "^    -      ''m  m 
^1  m 


yield  the  matrix  equation 


EG   g   =  X  (2-49) 

mn   n    m  ' 

n 

where 


,3. 


^mn  =   I   ^''^   \^L>    ^(L-i)  (2-50) 

K 

and 

\  =   i,^'r,'^\r)    Gjr)  (2-51) 


=  /d^r  p^/^r)  G  (r) 
•^         —   m  — 
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The  matrix  elements  G   may  be  evaluated  in  analytic  form  if  we  choose 
the  functions  G  (_r)  to  be  some  analytic  function  such  as  Hermi  te-Gauss  ian 
functions.   The  integrals  x  ,  however,  require  numerical  integration  for 
most  cases  where  only  the  charge  density  is  known  as  a  linear  combination 
of  analytic  functions.   This  procedure  does  allow  a  reduction  of  the 

total  number  of  numerical  integrations  required  compared  to  a  process 

2 
not  using  fitting  by  an  approximate  factor  of  N  /N  ,  where  N  is  the  number 

of  orbital  basis  functions  and  N   is  the  number  of  exchange  fitting  func- 
tions G  (r)  . 

m  — 

The  total  charge  density  P {r)    is  also  approximated  using  a  linear 
combination  of  analytic  functions  such  as  Hermi te-Gaussian  functions. 
The  rationale  for  fitting  the  charge  density  as  well  as  p^/^  is  that  the 
fitting   process  will  reduce  the  total  number  of  computationally  time 

consuming  integrals  necessary  for  computing  the  Coulomb  contributions  to 

2 
the  total  energy  and  the  secular  matrix  by  an  approximate  factor  of  N  /N  , 

where  N   is  the  number  of  charge  fitting  functions  F  (r).   Since  the  value 
C  3         3  ^\_/ 

of  N   is  typically  the  same  order  of  magnitude  as  N,  approximating  the 
density  may  produce  a  substantial  saving  in  computational  time. 

Before  explaining  the  procedure  for  choosing  the  coefficients  f  , 
let  us  digress  for  a  moment  to  the  problem  of  approximating  charge  den- 
sities in  the  general  case  using  the  methods  introduced  for  molecular 
claculations  by  Dunlap  (Dunlap,  Connolly,  and  Sabin,  1979a)  and  by 
Mintmire  (1979)-   Consider  a  prototypical  molecular  charge  density  o {r) 
such  that  0  decreases  more  rapidly  than  r    in  the  limit  of  large  r. 


Then  define  the  electrostatic  interaction  U  of  a  charge  density  with  it- 
2 


self  to  be  U  =  -T-  [p|p],  where 
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7  7  Pi  y  Po^L') 

[p,1p2]  =  /d\  jd\'         I   .  ^f| (2-52) 


If  the  specified  charge  density  p  is  approximated  by  another  density  p 
which  is  dependent  upon  several  parameters  which  may  be  varied  independ- 
ently, and  if 

p(£)  =  p(£)  +  Ap(£)  (2-53) 

then 

[pIp]  =  2[pI5]  -  [p|p]  +  [ApJAp]  (2-5A) 

Let  us  define  an  approximate  electrostatic  interaction  energy  U  in  terms 
of  the  approximate  density  p  and  the  exact  charge  density  p: 

U  =  [pIp]  -  J   [p|p]  (2-55) 

Since  the  difference  AU  =  U  -  U  =  -  [Ap|Ap]  is  greater  than  or  equal  to 


zero  for  any  choice  of  p,  a  reasonable  criterion  for  choosing  p  is  the 
minimization  of  AU.   This  constraint  will  be  satisfied  if  the  variation 
of  AU  with  respect  to  allowed  variations  of  p  is  zero: 

6AU  =  6U  =  0  (2-56) 

leading  to  the  relationship  expressed  in  terms  of  5p: 

<5U  =    [p|6p]    -    [pl6p]   =   0  (2-57) 

We  see  that  if  a  variation  of  the  form  6p  =  p  dc  is  possible  within 
the  range  of  variations  allowed  to  p  and  de  arbitrary,  then  the  additiona: 
relationship 
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{[pIp]  -  [p|p]}  de  =  0  (2-58) 
is  valid.   This  equation  implies  the  equality 

[p|p]  =  [pIp]  (2-59) 

which  shows  that  for  the  choice  of  p  satisfying  Equation  2-57,  then 

ij  =  J   [p|p]  =  J   [p|p]  (2-60) 

For  p  expressed  as  a  linear  combination  of  fitting  functions  F  (r) : 

m  — 

3(l)  =  ^  f^  ^^iL>  (2-61) 

m 

we  may  restate  Equation  2-57  in  terms  of  the  quantities  f  and  F  (r): 

m  m  — 


9p    ,         r~|9p 

af-J  "  fplgf 

mm  m 


6U  =  E    {[p||f-]    -    [p||f-]}    6f^  (2-62) 


=   E    {[p|F   ]    -   E   f      [F    |F    ]}    6f     =0 
m  n        n      m  m 

m  n 

Since  the  variations  6f  are  arbitrary  and  independent,  each  term  in  the 

m 

summation  over  the  index  m  must  equal  zero.   Rearranging  terms  yields: 

n 

or 


E  B   f  =  a  (2-64) 

mn  n    n 
n 

where 


Bn,n  =  ^^J^J  (2-65) 

mn     m '  n 


a  =  [F  |p] 
m     m' 
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With  the  proper  choice  of  fitting  functions  F^^  and  form  of  the  charge 
density  p,  these  integrals  may  be  evaluated  analytically. 

Extension  of  these  methods  to  extended  systems  requires  a  few  mod- 
ifications.  The  fitting  functions  are  now  the  periodic  fitting  functions 

F^(£)  formed  from  the  periodic  sum  of  the  localized  functions  F  (r)  as 

m  — 

stated  in  Equation  1-kl .      The  charge  interaction  expression  [pjp^]  must 
also  be  redefined  in  order  to  avoid  dealing  with  infinite  terms.   The 
most  convenient  definition  of  this  term,  if  p,  and  p^  are  periodic  extended 
charge  densities,  is  the  electrostatic  interaction  per  unit  cell: 

[pjp^]  ^V^r/d^r-   \   .  ^,,  (2.66) 


The  charge  interaction  term  [p^Ip^]  is  still  non-negative  for  p^  =  p 
and  pj  having  no  singularities.   In  addition,  this  expression  is  still 
symmetric  with  respect  to  the  order  of  p,  and  p^  (both  of  which  are 


periodic),  since 

>,|P2]  =  /^a  r  jd-r'   i   .  ^.T (2-67) 


r   I   T    f  3   ^  3    P]  (Jl)  P?(l') 
[pjpj  =  Ld\  /d\'   I   _  ,.?- 


3   ,   ,    0,(£)  0  (r'  -  R) 
=  z//r/^d3r'-I- ^-    - 


R 


r  -  r'  +  R 


Pi^L  +  R)  Oo(r') 


3   r       ■\  ■l^.L^H^'^oVr 


r  +  R  -  r' 


,3,3  P2(L')  Pi^"-) 

=  V  '■'//'■  ^ '—     =    fP^lo,] 

r  -  r ' 


This  expression  for  the  charge  interaction  is  also  convenient  because 
It  satisfies  two  conditions: 
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0   for  ^^(r)    a    periodic  function  formed  from  the  periodic  sum  of 
of  localized  functions  F  (r)  as  in  Equation  Z-^*?,  then 


[F  If  ]  =  1  j6\  j6\'    '"^-'  ■"^-  "' 


F_(r)  F  (r  -  R) 

(2-68) 


m'  n    ^   '  '  IL  -  L  I 

which  is  analytically  integrable  in  closed  form  for  F  (r)  being 

m  — 

a  Hermi te-Gaussian  function; 
2)   [pjp^]  is  finite  if  and  only  if  one  or  both  of  the  charge 

distributions  p.  and  p-  satisfy  the  charge  neutrality  condition: 

\^^r   p(£)  =0  (2-69) 

The  first  condition  is  evident  upon  inspection.   The  second  condition 
follows  from  consideration  of  a  limiting  procedure  which  begins  with  a 
finite  assembly  of  identical  unit  cells.   The  initial  assembly  is  con- 
structed with  no  embedded  voids.   The  limiting  process  then  consists  of 
arbitrary  enlargement  of  the  finite  assembly,  with  the  infinitely  periodic 
(in  two-dimensions)  slab  the  ultimate  result.   Any  such  finite  assembly 
has  a  periodic  charge  density  which  satisfies  the  explicit  constraint  on 
asymptotic  behavior  required  by  the  definition  of  the  molecular  charge 
density  interaction  [p.|p-]  of  Equation  2-52.   It  then  follows  from 
Equation  2-52  that  the  charge  interaction  for  such  a  finite  assembly  with 
the  charge  density  belonging  to  a  s  ingle  unit  cell  is 

-     -   p,  (0  p  (r ') 
R_  II  "  L  ^  ^1 

where  the  primed  summation  indicates  a  finite  number  of  terms.   In  order 
to  carry  through  a  physically  reasonable  limiting  process  the  single 
unit  cell  is  conceived  as  being  the  most  nearly  central  unit  cell  of  the 
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finite  assembly.   Similarly  the  origin  of  coordinates  is  associated,  for 
convenience,  with  the  central  unit  cell.   Then  for  each  term  of  Equation 
2-70  save  for  the  one  corresponding  to  R^  =  £  the  contribution  may  be 
exactly  as  a  multipole  expansion  (see  Appendix  B  for  a  more  detailed 
discussion).   if  the  charge  neutrality  condition  is  not  satisfied  then 
the  leading  term  in  the  multipole  expansion  will  be  a  monopol e-monopole 
interaction  term  proportional  to  R~  .   For  a  two-dimensional  lattice  of 
vectors  R^  an  infinite  sum  of  monopole-monopol e  interactions  is  not  con- 
vergent.  Therefore  in  the  limit  of  an  arbitrarily  extended  slab,  the 
charge  neutrality  condition  is  necessary  if  the  charge  interaction  is  to 
converge  to  finite  limit.   Similarly  if  the  charge  neutrality  condition 
is  satisfied,  the  leading  term  in  the  multipole  expansion  for  given  R 
is  a  monopole-dipole  term  proportional  to  r"^.   It  is  shown  in  Appendix 
B,  however,  that  this  multipole  term  has  zero  total  contribution  to  the 
charge  interaction  term  for  two-dimensional  lattices  with  inversion  sym- 
metry.  Thus  the  leading  order  term  with  a  non-zero  contribution  to  the 
charge  interaction  term  will  be  a  term  proportional  to  R"^,  which  leads 
to  a  convergent  sum  in  Equation  2-70  in  the  infinite  limit  under  consid- 
eration. 

Thus  the  second  condition  leads  immediately  to  the  observation  that 
Equations  2-63  through  2-65  as  derived  are  not  usable  in  a  direct  manner, 
since  the  density  involved  is  purely  electronic  and  manifestly  not  neutral 
However,  Equation  2-65  is  valid  for  the  extended  system  so  long  as  charge 
densities  neutral  by  cells  are  employed.   For  any  finite  system  Equation 
2-63  may  be  rewritten  trivially  as 

^'>    -  ^  ^n  ^1  =  0  (2-71) 

n 
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irrespective  of  charge  neutrality.   In  the  limit  described  above,  however, 
the  charge  interaction  expression  becomes  that  of  Equation  2-67  and  the 
necessary  and  sufficient  conditions  on  charge  neutrality  come  into  play 

with  the  result  that  if  at  least  one  of  the  F   is  not  charge  neutral, 

m 

then  the  right-hand  factor  in  Equation  2-71  must  be  charge  neutral, 
where 

/.d^r  p(r)  =  E  f   f  d^r  F  (r)  =  S  f  jd^r   F  (r)  (2-72) 

^  u  —  n'i2     n—       n^     n  — 


n 


or  expressed  differently 


Z  f  n  =  1  (2-73) 

m  m 
n 

where 

n  =  {jd^r    F  (r)}  /  {/^d^r  p(r)}  (2-7^) 

m         m  —      '  iZ      — 

We  can  produce  a  matrix  equation  equivalent  to  Equation  2-6^4  if  we  intro- 
duce a  nuclear  lattice  charge  P^.  to  each  electronic  density,  thereby 
giving  the  requisite  cell-by-cell  neutrality,  to  the  interaction  term  in 
Equation  Z-Gk   and  then  subtract  out  terms  contained  in  Equation  2-71  and 
analyze  the  identity 

[F     -   n  P,,|p   -   P^]    -    [F^  -   n  p^|e   f    (F     -   n^P^)]  (2-75) 

m  mN  N  m  mN  nn  nN 

n 

=    [F    Ip    -    E    f    F    ]    -    {1    -    E    f   n    )[F      -   n   o,,|p.,]    -    n^lp^lp    "    E    f    F    ] 
■•   m'  nn  nnm  mN'N  mN'  nn 

n  n  n 

The  first  term  of  the  right  hand  side  of  Equation  2-75  is  zero  according 

to  Equation  2-71.   Since  Equation  2-73  states  that  (1  -  E  f  n  }  =  0  and 

n   ''  ^ 

in  addition  we  see  that  [F   -  n  P.,|pn]  is  finite,  then  the  second  term  of 

m    mN'N 
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the  right-hand  side  of  Equation  2-75  must  equal  zero.   The  term 

[p  |p  -  X  f  F  ]  must  also  be  finite  (again  by  charge  neutrality)  for 
N      f^   n  n  =>  I  ^  II 

thin  films;  therefore  Equation  2-75  reduces  to: 


{E  f  [F 
n  m 

n 


m  N '  n 


n  N 


-  [F 


m 


m  N ' 


N     m 


(2-76) 


where 


X  =  [pjp  -  E  f^Fj 
n 


(2-77) 


is  an  undetermined  constant  independent  of  the  index  m  in  Equation  2-76, 
Restating  Equation  2-76  in  matrix  notation  yields: 


Z  B   f 

mn   n 
n 


a  +  Xn 
m     m 


(2-78) 


where 


B   =  f"^  -  ^  PmI"^  -  n  p.,] 
mn     m    m  N'  n    n  N-" 


(2-79) 


^r^   ~    f'^     "   ^   PmIp   "   PkJ 

m     m    m  N '     N 


(2-80) 


This  set  of  equations  is  similar  to  the  equations  resulting  from  a 
least-squares  fitting  procedure  which  includes  a  linear  constraint  equation. 
Thus  we  may  solve  for  A  by  rearranging  Equation  2-78: 


f   -  E  {B"b   a  =  A  E{B"b   n 
m   ^      mn   n    „     mn   n 
n  n 


(2-81) 


E  n^  f^  -  I  E  n  {B"h   a  =  x  e  E  n  {b"'}   n 

mm        m      mn   n        m      mn   n 
m         m  n  m  n 


resulting  in  the  value  of 


(2-82) 


A  =  {1  -  E  z  n„  {B"h   a  }  /  {E  E  n  {b"S   n  } 
m      mn   n         m      mn   n 
m  n  m  n 


(2-83) 


3h 


The  solution  for  the  coefficients  f   is  then  given  by 

f  =  r  {B"h   -(a  +  An  )  (2-8A) 

m         mn   n     n  ^  ' 

n 

We  see  that  the  approximate  electrostatic  interaction  energy  defined  in 

Equation  2-55  now  has  the  form: 

U=2:f  a  -  ]r     ZZf   {b"^}   f  (2-85) 

n  n   2       n      mn   n 
n  m  n 

=  [p!p]  -  J   [5|p]  -  [p|p^]  +  J   [pjpj  (2-86) 

which  is  the  electrostatic  energy  of  the  combined  electronic  and  nuclear 
charges  with  only  the  electron  repulsion  term  being  replaced  by  an  approx- 
imate interaction. 


2-5.  Self-Consistent  Solution  of  the  Secular  Equation 

The  orbital  coefficients  c  •(]<//)  are  obtained  by  solving  the  secular 

equation  in  Equation  2-'*l.   Since  the  matrix  elements  H   (k,,)  are  de- 

pendent  upon  the  coefficients  c  •(]<.//)  for  which  we  seek  a  solution, 

Equation  2-^1  is  solved  using  the  familiar  iterative  procedure  known  as 

the  self-consistent  field  (SCF)  method,  described  by  Slater  (1963)  for 

atomic  and  molecular  systems.   In  the  SCF  method  one  starts  with  an 

initial  choice  of  c  .(k,,)  to  generate  the  matrix  elements  H   (k,,).   The 

ni  — //  mn  — // 

secular  equation  is  then  solved  to  yield  a  new  set  of  values  for  the 

coefficients  c  .(k,,)  which  is  used  to  begin  a  new  iteration.   This  process 
n  I  — // 

is  repeated  until  some  predetermined  criterion  of  convergence  for  the 
iterative  process  is  satisfied. 
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Due  to  the  introduction  of  the  fitting  procedures,  the  secular 
equation  defined  in  Equations  2-41  through  2-43  will  not  use  the  effective 
Hamiltonian  as  it  stands  in  Equation  2-34.   Instead  the  effective  Ham- 
iltonian  will  be  of  the  form  (for  the  spin  up  effective  Hamiltonian): 

H+.r(r)  =  -i  V^  +  E  f   /d^'  -r-^ r  -  Z    E  n f -J- 

eff  —     2         m-'      r-r    „    r-R-R 
m  '      '    R  a  '—   — a   — ' 


3il/3 

m  m  — 


-  3ct[7;7r''  ^   g„  G„(r)  (2-87) 


m 

Thus  given  an  initial  choice  of  c  •(]<//)  and  n.(j<,,),  the  logical 
order  of  each  iteration  will  be: 

1)  Compute  the  fitting  coefficients  f  and  g  ; 

m      m 

2)  Use  the  f i tt ing  coefficients  to  generate  the  secular  matrix 
defined  in  Equation  2-87; 

3)  Solve  the  secular  Equation  2-41  to  yield  a  new  set  of  coefficients 
c  .  (J<  / /)  and  occupation  numbers  n.(k^,,). 

Since  this  procedure  differs  from  those  of  typical  LCAO  calculations  In 

the  use  of  a  fitting  formalism,  we  present  briefly  a  discussion  of  the 

necessary  equations  for  the  above-mentioned  procedure  for  each  iteration, 

beginning  with  an  Initial  set  of  coefficients  c  •(]<//)  and  occupation 

numbers  n  .  (k  ,  ,)  . 
I  — // 

First  the  charge  density  p  of  Equation  2-35  is  used  with  the  fitting 
scheme  of  Equations  2-78  through  2-84  to  compute  the  charge  fitting 

coefficients  f  .   The  B   matrix  elements  of  Equation  2-79  are  evaluated 

m        mn 

directly  while  the  values  for  a   are  constructed: 

m 

a  =  {E  E   n  (k   )  E  E  c''.\(k.,)  c,  .(k..)  V.,  (k,,)}  -  C     (2-88) 
m    ;  K    I  — //   .  ,   J  I  -//   ki  -//   jkm— //'     m 
_//         J  "^ 
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where 


Vjk^(Ji//)    =   ^    ^      exp[ik^/-R]    /d^r   U?(£+   R)    U|^(r_) 
R    R 


,        F    (r'    -   R')  2 

X   (/«'r'      J,.,,  I  -   n„  I    I,,    ,%     ■„.,)  (2-89) 

'—       —   '  a    '—         ^a        —   ' 


and 


a        1   _      r-   T- 1         a      b 


R_  a    '—       — a'  a    b    '—a       — b       — ' 

The  primed  summation  indicates  that  nuclear  sel f- i nteract ion  terms  for 

which  R^  is  the  zero  vector  and  indices  a  and  b  denote  identical  nuclear 

sites  are  to  be  omitted. 

Next  the  cube  root  of  the  density  is  fit  using  Equation  1-k'^ .      The 

G   matrix  elements  are  evaluated  in  the  form  given  in  Equation  2-50. 
mn  3         -1         ^ 

The  X   integrals  are  evaluated  using  a  numerical  integration  procedure. 

For  each  numerical  grid  point  r_.    the  value  of  the  density  P  ( r_. )  is  generated 

from  the  precalculated  values  of  B1och  functions  ii;.(r.;  k,,)  and  the 

J  —I   -// 

exchange  functions  G  (r.): 
^  m  — I 

p(r.)  =  E   E   n;,(k..)  ]  1   c.«(k,/)  ^.{r.\    k.,)  1^  (2-91) 

-//  -• 

which  is  then  used  to  compute  p^^(_r.).   The  numerical  cube  root  of  the 

density  is  used  then  v^ith  the  values  of  the  exchange  fitting  functions 

G  (r.)  to  evaluate  x   integral  numerically: 
m  — I  m     ^  ' 

x^  =  E  w.  p^^^(r.)  G^(r.)  (2-92) 

m    .   I      — I   m  — I 
I 

where  the  weights  w.  and  the  grid  points  r_.    are   appropriately  chosen 
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grid  points  for  a  three-dimensional  numerical  integration  over  the  unit 
eel  1  vol ume  fi. 

The  fitting  coefficients  f  and  q   are  then  used  to  evaluate  the 

m     -m 

elements  of  the  secular  matrix  defined  by  Equation  2-^3.   Using  the 
expression  for  the  effective  Hamiltonian  of  Equation  2-87,  the  matrix 
elements  for  the  spin  up  effective  Hamiltonian  have  the  form: 

H..(k/,)  =  T..(k,.)  +  If  v..  (k..)-3a[-4]^^^j:q  X..  (k,,)   (2-93) 
I J  — //     I  J  — //      m  ijm  — //   ^   i»7r      "m  ijm  — // 

where 

T..(k,,)  =  E  exp[ik.,-R]  /d^r  UV  (^  +  R)  {-y  V^}  IJ .  (O         (2-9'+) 

IJ//n  /'  '  ^         J 

and 

X..  (k  )  =  Z   exp[ik..-R]  y      /d^r  U^(r  +  R)  U.(r)  G  (r  -  R')  (2-95) 
ijm  -//    ^  -II   -f^i       "-   -   J-   m-   - 

This  form  of  the  secular  matrix  is  diagonal ized  and  new  occupation  numbers 

n.(j<.  ,)  are  computed,  using  the  eigenvalues  of  the  secular  equation  to 

find  the  lowest  eigenstates  of  the  effective  Hamiltonian  to  occupy.   At 

this  point  the  agreement  between  the  current  orbital  coefficients  c  .(k,,) 

ni  -// 

and  those  present  at  the  beginning  of  the  current  iteration  is  tested  for 
convergence  of  the  iterative  sequence.   If  the  convergence  criteria  are 
not  satisfied,  a  new  iteration  begins. 


CHAPTER  I  I  I 
COMPUTATIONAL  METHODS 


3~1.  Choice  of  Basis  Sets 

In  any  atomic,  molecular,  or  solid  state  calculation  of  electronic 
structure  using  a  linear  combination  of  a  finite  set  of  basis  functions 
as  approximate  orbitals,  tine  particular  choice  of  functions  is  a  critical 
part  of  the  computational  procedure.   The  first  decision  in  choosing 
basis  functions  that  must  be  made  is  the  type  of  functions  to  be  used. 
This  work  is  based  on  the  use  of  Hermi te-Gauss ian  functions,  which  are 
described  in  Appendix  A.   These  functions  are  specific  linear  combinations 
of  the  more  familiar  Cartesian  Gaussians  of  the  form  x  y  z  exp(-ar  ), 
where  1,  m,  and  n  are   nonnegative  integers  and  the  exponent  coefficient 
a  (commonly  denoted  the  "orbital  exponent")  is  constant  for  all  Cartesian 

Gaussians  which  are  components  of  the  Hermi te-Gauss ian .   That  is,  all 

2 
Hermi te-Gauss ians  may  be  expressed  as  P (x,y , z)exp(-ar  ),  where  P(x,y,z) 

is  a  polynomial  of  finite  order  of  the  coordinate  variables  x,  y,  and  z. 
The  total  number  of  Hermi te-Gauss ian  functions  in  the  basis  and  the 
choice  of  the  orbital  exponents  a.  is  also  a  critical  choice.   Fortu- 
nately there  is  an  extensive  literature  of  results  of  atomic  and  molec- 
ular calculations  using  Cartesian  Gaussian  functions  that  may  be  consulted 
in  making  these  decisions  so  that  Cartesian  basis  sets  may  be  converted 
to  Hermi te-Gauss ian  basis  sets.   The  conversion  for  most  systems  from 
Cartesian  Gaussian  basis  sets  to  Hermi te-Gauss ian  basis  sets  is  usually 
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trivial,  since  s  and  p  type  Cartesian  Gaussian  functions  correspond  to 
equivalent  Hermi te-Gaussian  basis  functions.   Cartesian  d  type  functions 
cannot  be  expressed  as  a  single  Hermi te-Gauss Ian  function;  instead  they 
must  be  expressed  as  a  sum  of  two  Hermi te-Gauss ian  functions.   This  is 
not  a  serious  problem,  since  we  may  include  both  Hermi te-Gauss ian  func- 
tions and  contract  them  to  yield  the  Cartesian  Gaussian.   Thus  we  may 
use  Cartesian  Gaussian  basis  sets  expressed  as  contracted  Hermite- 
Gaussians  if  necessary.   The  charge  fitting  and  exchange  fitting  basis 
sets  are  then  chosen  using  the  method  described  by  Dunlap,  Connolly,  and 
Sabin  (1979a)  for  atomic  and  molecular  calculations. 

3-2.  Overlap  and  Kinetic  Energy  Integrals 

The  overlap  matrix  elements  S.  (_k^^)  between  normalized  Hermite- 
Gaussian  basis  functions  will  be  of  the  form: 

K 

/d^  f(fii,a.,A.-R;  r)  f{n.,a.,^.;    r)  (3-I) 

where  N.  and  N  are  normalization  constants  for  the  localized  Hermlte- 
Gaussians.   The  kinetic  energy  matrix  elements  T.  .  (J<  , -)  are   evaluated  at 
the  same  time  as  the  overlap  matrix  elements  since  they  are  of  the  same 
general  form  as  the  overlap  matrix  elements  such  that: 

T..(k//)  =  (N.N.)"^/2^exp[ik//-R] 

R 

1    ^     3 

2   Z  /d\  f(n  +2G  a  ,A  -R;  r)    f(n.,a.,A.;  r)  (3-2) 

S=l  ill  J   J   J 


ko 


where  the  integer  sets  u   refer  to  the  sets  (1,0,0),  (0,1,0),  and  (0,0,1) 

for  values  of  s  equal  to  1,  2,  and  3  respectively.   For  values  of  a  and 

2 

b  such  that  (a+b)R    is  greater  than  1  (the  value  of  1  is  rather  arbi- 

max 

trary  and  may  be  adjusted  in  order  to  speed  computational  time  as  more 
computational  benchmarks  are  accumulated)  these  integrals  will  be  eval- 
uated using  direct  lattice  sums  of  integrals  of  the  form  given  in  Equa- 

2 
tion  A-5  in  Appendix  A.   For  values  of  (a+b)R    (where  R    is  the 

max         max 

greater  of  R^.  and  R^^'  the  primitive  lattice  vectors)  less  than  or  equal 
to  1,  the  overlap  matrix  elements  and  kinetic  energy  matrix  elements  are 
evaluated  using  reciprocal  lattice  space  techniques  described  in  Appendix 
C  using  expansions  of  the  form  given  by  Equation  C-8. 

Lowdin  orthonormal izat ion  (Lowdin,  1956)  is  used  to  transform  the 
matrix  elements  of  the  secular  matrix  to  an  orthonormal  basis  so  that  the 
S  matrix  in  the  secular  equation  will  be  the  unit  matrix.   First  the 
unitary  matrix  U  is  found  which  diagonal izes  the  overlap  matrix: 

y^  s  y  =  d  (3-3) 

The  overlap  matrix  may  then  be  transformed  to  the  unit  matrix  by  use  of 
the  matrix  T: 

r   S  T  =  1  (3-4) 


where 

T  =  U  d"^/2  (3_3) 


The  secular  equation  of  Equation  2-^41  then  transforms  according  to 

t'  H  T  (t"^  C)  =  (t"^  C)  e  (3-6) 
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Approximate  1  Inear  dependencfes  (which  arise  from  the  limits  to  precision 
inherent  in  any  computational  procedure)  in  the  basis  may  be  handled  by 
excluding  eigenvectors  of  the  overlap  matrix  v-^hich  have  eigenvalues  less 
than  some  specified  tolerance  (for  the  purposes  of  this  work  the  tolerance 
was  chosen  to  be  1 0  ). 

3~3-  Fitting  Function  Integrals 


The  integrals  B..,  C.,  V..  (k.,),  and  X..  (k,,)  currently  are  eval- 
1 J    J    1  jm  — //        I  jm  — // 

uated  using  direct  lattice  expansion  techniques  as  discussed  in  Appendix 
A.  The  matrix  elements  B..,  defined  in  Equation  2-79,  may  be  reduced  to 
the  form: 

n.       5/2        ^   ^    a. a. 

B..  =  Z{(-1)  J  ^ — -  R(n.+n.,    '  ■>  ,  A.-R-A.) 

'J   R       a.a.(a.-Ha.)^/2     .   j   ^  .^^  _   -.  - -j 

'  J   '   J  'J 

n.  n. 

-  Z    2ttI    [   —  R(n.,a.  ,A.+R-R  )  +  — L  R(n .  ,a .  ,A. -R-R  )] 
^    a  a     J   J  -J a    a .     i   i  '-i a 

•^  a  b   '—a -h ' 


using  Equations  A-^0  and  A-kS.      The  C.  integrals,  defined  in  Equation 
2-90,  may  be  reduced  to 

C   .  r  (  z  !:!l  R(S   a   A  -R-R^)  -  n   I  E'  j^^   )  (3-8) 

-•   Raj      jjj         -"ab'— a -b' 


Equations  A-'+l  and  A-A7  may  be  used  to  convert  the  matrix  elements 

v.  .  (k  ,,)  to  the  form: 
ijm  — // 


k2 


V::^('i//)  =  (N.N.)'^^^  I    exp[ik   -R]  I   g  (K,  n .  ,n  .  ,a  .  ,a  . ) 


a. a.  n     «  5/2 

f(K.— ^-^,A.-R-A.;  0)  1    {(-1)  "^    ^^ 


a.+a.  — 1 1   —  „i        /     \   /        \l/2 

■   j       -•     R'       (a.+a.)a  (a.+a.+a  ) 
—  I   J   m   1   J   m 


(a.+a.)a  a  .  (A.-R)+a .A . 

R(n.+n.+n  -K,    J       '  "^ ,  '  "'  ~  ^~^    -   A  -R') 

I   J   m    a.+a.+a  a.+a.        -m  — 

I   J   m  I   J 


2TTn  a.  (A. -R^)+a.A. 

■1    1  ^  R(n.+n.-K,  a.+a.,   — — ^ ^-L  -  r  -r')} 

a  a.+a.     i   j  —  i   j    a.+a.  -a  — 

a     I   J       ■'  ■'  I   J 


(3-9) 

For  the  direct  lattice  expansions  displayed  in  Equations  3-7  and  3-8  and 
the  inner  lattice  sum  over  R^'  in  Equation  3-9,  multipole  expansion  tech- 
niques as  described  in  Appendix  B  are  used  to  increase  the  computational 
speed  of  evaluating  these  integrals. 

The  exchange  fitting  integrals  X..  (k,,),  defined  in  Equation  2-95, 

ijm— // 

are  also  evaluated  in  a  double  lattice  sum  having  the  form: 

X..  (k,,)  =  (N.N.)"^'^^E  exp[ik,/R]  I  g(K,n.  ,n  .  ,a  .  ,a  . ) 
ijm  -//      I  J      R       -^^  ~  K      '   J   '   J 


I   J        ■'  R.    '   J   "^ 


(a.+a.)a   a . (A. -R) -a .A . 

f(n.+n.+n  -K,    '   {  ^  ^^^ ^  -  A  -R';  0) 

I   J   m    a.+a.+a  a.+a.  — m  — 

-•         I   J   m       I   J 


(3-10) 
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All  of  the  above  Integrals  could  also  be  evaluated  using  reciprocal 

lattice  expansions  as  discussed  in  Appendix  C.   Since  V..  (k,,)  and 

I jm  ^/ 

X..  (k,,)  each  contain  two  lattice  sums  as  expressed  in  Equations  3"9 
1  jm  — // 

and  3"10,  we  may  choose  to  convert  one  or  both  of  the  direct  lattice 

summations  to  reciprocal  lattice  expansions.   An  approach  using  double 

reciprocal  lattice  expansions  has  been  tested  for  the  case  of  all  sphei — 

ically  symmetric  Gaussian  functions  for  the  V. .  (k,,)  and  X..  (k,,) 

ijm  — //       ijm  — // 

integrals  (which  are  the  computationally  most  time  consuming  integrals 
to  evaluate).   The  results  from  these  tests  indicate  that  the  large  number 
of  sine  and  cosine  function  evaluations  required  make  the  double  lattice 
reciprocal  lattice  expansion  algorithm  at  best  equal  in  computational 
time  to  the  direct  lattice  expansion  algorithm  for  integrals  involving 
diffuse  Gaussian  functions,  which  is  the  case  requiring  the  most  computer 
time  using  the  direct  lattice  expansion  algorithm.   For  Gaussian  functions 
which  were  not  extremely  diffuse,  the  double  reciprocal  lattice  expansion 
algorithms  required  considerably  more  computational  time  than  algorithms 
using  the  double  direct  lattice  expansions. 

Another  approach  would  be  to  convert  only  the  inner  lattice  sums 
over  R^'  in  Equations  3~9  and  3"10  to  reciprocal  lattice  expansions.  The 
formulas  necessary  for  incorporating  this  approach  into  computational 
algorithms  are  derived  in  Appendix  C.  This  approach  has  not  yet  been 
incorporated  into  the  computer  code,  leaving  open  one  avenue  for  increasing 
the  computational  speed  of  the  computer  code  for  calculating  the  thin 
layer  electronic  structure. 
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3-A.  Numerical  Integration  Techniques 


In  order  to  fit  the  exchange  potential  [or  more  precisely,  the 

1/3 
function  P    (r)]    to  a  linear  combination  of  fitting  functions  G  (r) 

m  — 

we  must  numerically  integrate  the  expression: 

x,  =  /,  d^P^/^Cr)  G^(^)  (3-n) 

This  leads  to  the  approximate  form  of  x  : 

m 

\~-   ^-ip'^'(Li)  Gjr.)  (3-12) 

I 

The  integration  points  £.  and  weights  w. ,  which  are  chosen  for  numerical 

integration  over  the  para  1 lelopiped  primitive  unit  cell  defined  in 

Section  2-1,  are  discussed  in  Appendix  D. 

The  values  of  the  primitive  orbitals  ■•>! .  {r^- ;    k,,)    and  the  values  of 

the  periodic  fitting  functions  G  (r.)  must  be  evaluated  at  each  inteqra- 

m  — I  ^ 

tion  point.   For  a  primitive  Bloch  orbital  ^{r_.;    U^..)    constructed  from 
localized  Hermi te-Gauss ian  functions: 

■^(L;;  i//)  =  n'^^^  T   exp[ik^^-R]  f(n,a,A+R^;  £. )  (3-I3) 

K 

we  may  evaluate  i|'(_r.  ;  k^..)    either  by  using  the  above  direct  lattice 

expansion  or  by  using  the  reciprocal  lattice  expansion: 

"'('"!'  "<//)  =  h   n"''''  a"^    H   [a'^'-'Cz-A  )] 
— i      —/ /  au  n,         2 

n^+n,  n^  n, 

^    {-<)  ^    (K  +k.,  )    (K  +k,,  )  ^ 

1^  x  -//x      y  -//y 

exp[-(K+j^//)  /ka]    exp[i(K+k//)-(£.-A)]  (3-1^4) 
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discussed  in  Appendix  C.   Evaluation  of  the  periodic  fitting  function 

G  (r.),  defined  in  Equation  2-46,  is  equivalent  in  form  to  that  of  a 
m  — I 

Bloch  function  with  a  wave  vector  j<,,  equal  to  the  zero  vector.   Thus 
we  may  evaluate  the  fitting  function  with  either  a  direct  lattice  expan- 
sion or  a  reciprocal  lattice  expansion: 


G(_r  )  =  I   f(n,a,A  £  )  (3-15) 

=  -TT  a     H   [a  '  (z-A  )]  E  (-r) 
af2        n3         z   ^ 

"l      "2  2 

K        K    '   exp[-K  Aa]    exp[  i  K- (j^.-A)  ] 
X       y  I 


The  matrix  elements  G   ,  defined  in  Equation  2-50,  may  be  evaluated 

mn 

either  analytically  or  numerically.   We  have  chosen  to  evaluate  G 

mn 

numerically  using  the  same  integration  grid  as  used  for  the  evaluation 

of  X  .   This  particular  choice  of  the  same  numerical  integration  grid  is 
m  3       3 

advantageous  since  the  use  of  G    in  Equation  2-^49  evaluated  as 

mn     ^ 

G   =  Z  w.  G  (r.)  G  (r.)  (3-16) 

mn    .   I   m  — I    n  — I 
I 

is  equivalent  to  a  least-squares  fitting  procedure  over  a  weighted  discrete 

set  of  points.   That  is  to  say,  the  use  of  G   as  it  is  expressed  in 

mn 

Equation  3-16  in  the  solution  of  Equation  2-^9  yields  a  set  of  coeffi- 
cients for  the  exchange  fitting  functions  q  which  minimize  the  quantity: 

m 

A  =  E  w.  (p^^^j:;)  -  E  g^  G^^,))  (3-17) 

i  m 

V/e  find  that  this  particular  means  of  evaluating  G   yields  a  method 

mn 

which  requires  fewer  integration  points  (and  thus  less  computational 
time)  for  the  same  precision  in  the  exchange  fitting  procedure  than  if 
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we  evaluate  G   analytically.   This  is  not  difficult  to  understand  when 

one  considers  that  the  process  which  uses  the  numerically  evaluated 

values  of  G   will  solve  the  least-squares  fitting  procedure  over  a  dis- 
mn  ->  3  r 

Crete  set  of  points  with  an  error  of  the  order  of  the  precision  of  the 

evaluated  values  of  ii.{r.;    k , ,)    and  G  (r.),  which  are  typically  of  the 

J  —I   — //       m  — I 
_r      _g 

order  10   to  10   .   The  use  of  an  analytically  evaluated  G   matrix 

mn 

yields  a  process  which  will  have  an  error  of  the  order  of  the  precision 

of  the  numerical  integration,  which  is  typically  of  the  order  10   to 

-3 
10   .   As  long  as  the  number  of  grid  points  exceeds  the  number  of  fitting 

functions  G  (r)  such  that  the  discrete  least-squares  fitting  procedure 
m  —  ^ 

is  not  capable  of  yielding  solutions  with  severe  oscillations  between 
the  discrete  set  of  grid  points,  it  appears  that  the  use  of  a  numerically 

integrated  value  of  G   is  preferable. 

^  mn    ^ 


CHAPTER  IV 
RESULTS 


^-1.  Atomic  Hydrogen  Monolayer 

The  atomic  hydrogen  monolayer  was  chosen  as  the  initial  test  case 
for  the  approach  outlined  in  Chapters  2  and  3  for  calculating  the  elec- 
tronic structure  of  thin  films  because  atomic  hydrogen  possesses  only 
one  electron  and  may  be  treated  adequately  using  only  a  small  number  of 
basis  functions.   Calculations  have  been  performed  on  two  different  lat- 
tice structures:   a  square  lattice  and  a  hexagonal  lattice.   Preliminary 
calculations  were  performed  (Mintmire  and  Sabin,  igSOb)  using  a  minimal 
basis  set  consisting  of  three  s-type  Gaussian  functions  taken  from  van 
Diujneveldt  (1971)  for  the  orbital  basis.   The  basis  set  used  for  fitting 
the  charge  density  contains  three  s-type  Gaussian  functions  with  exponents 
equal  to  twice  the  exponents  used  in  the  orbital  basis.   The  basis  set 
used  for  fitting  the  exchange  potential  similarly  contained  three  s-type 
Gaussian  functions  with  exponents  equal  to  one-third  the  exponents  used 
in  the  charge  fitting  basis.   This  set  will  henceforth  be  referred  to  as 
Basis  1.   Table  4-1  presents  the  orbital  exponents  used  in  the  orbital 
basis  set  as  well  as  the  orbital  exponents  for  the  Gaussian  functions  used 
to  fit  the  charge  density  and  the  exchange  potential. 

The  value  of  2/3  was  used  for  a    in  all  calculations  in  this  work, 
except  where  explicitly  noted.   Other  values  of  a  could  of  course  be  used, 
such  as  the  a^^  or  a^^   values  for  the  atomic  systems  given  by  Schwarz  (1972) 
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Table  A-1.  Basis  sets  for  the  atomic  hydrogen  monolayer, 


Type  of  basis 
Orbital 


Charge  fitting 


Exchange  fitting 


Funct  ion 
Type 

s 
s 

s 

X 

y 

z 

s 

s 

s 

2^  2 
X  +y 

s 

s 

s 

1^   2 
X  +y 


Exponents 
Basis  1  Basis  2 


0.151398 
0.68IA40 
4.50180 


0.302796 
1.362888 
9.00360 


0. 100932 

0.'45'4296 
3.00120 


0.151398 

0.68^140 

4.50180 

1.000 

1  .000 

1.000 

0.302796 
1.362888 
9.00360 

2.000 

0. 100932 

0.454296 
3.00120 

1.666667 


The  nomenclature  for  the  function  type  is  of  the  form  x  V  ^z  ^  for 
non-spherical ly  symmetric  Hermi te-Gaussians.   The  superscripts  nj,  n2, 
and  n3  refer  to  the  set  of  integers  r\   =  (nj,  n2,  n3)  which  define  the 
Hermi te-Gauss ian  function  and  do  not  represent  the  Cartesian  Gaussian 
function  one  would  normally  associate  with  this  nomenclature.   Expres- 
sions of  the  form  "x^+y'^"  are  to  be  interpreted  as  the  sum  of  two 
Hermi  te-Gauss  ian  functions  of  type  x^  and  y-^. 


^3 

The  formalism  described  in  this  work  assumes  that  the  value  of  a  is  con- 
stant over  all  regions  of  space  as  opposed  to  the  muffin-tin  multiple- 
scattering  Xa  method  described  by  Johnson  (1973),  where  a   may  have  dif- 
ferent values  in  different  regions  of  space.   Since  one  of  the  ultimate 
goals  of  this  method  is  to  describe  the  electronic  structure  of  systems 
containing  more  than  one  type  of  atomic  center,  we  feel  that  it  is  simpler 
at  this  time  to  avoid  the  complicated  question  of  what  value  of  a  is 
appropriate  for  a  particular  system  and  use  the  value  of  a  resulting  from 
the  original  derivation  of  the  local  density  functional  formalism  approach 
of  Kohn  and  Sham  (1965)  and  Gasper  (195^). 

The  values  of  the  v;ave  vector  k_,,  were  chosen  such  that  all  orbitals 
were  periodic  over  a  translation  equal  to  four  times  that  of  any  primitive 
lattice  translation,  which  is  equivalent  to  a  Born-von  Karman   region 
equal  to  a  four  by  four  array  of  contiguous  primitive  unit  cells.   Thus 
this  choice  yields  sixteen  evenly  spaced  non-equivalent  (by  translation 
symmetry)  wave  vectors  in  the  central  Brillouin  zone.   For  the  square 
lattice,  point  group  symmetry  reduces  the  number  of  non-equivalent  wave 
vectors  to  seven.   The  Brillouin  zone  of  the  square  lattice  is  depicted 
in  Figure  '4-1,  with  the  irreducible  region  of  the  Brillouin  zone  defined 
by  the  triangle  (rXM) .   Figure  k-2    illustrates  the  location  of  the  seven 
non-equivalent  wave  vectors  used  in  the  calculations.   The  hexagonal  lat- 
tice has  only  four  non-equivalent  wave  vectors  in  the  Brillouin  zone. 
The  hexagonal  Brillouin  zone  and  the  four  non-equivalent  wave  vectors  are 
Illustrated  in  Figures  '♦-B  and  k-h.      The  cohesive  energies  per  particle 
were  calculated  by  taking  the  difference  of  the  total  energy  at  a  given 
lattice  constant  and  the  energy  of  atomic  hydrogen  calculated  using  the 
molecular  LCAO-Xa  program  (Mintmire,  1979).   These  results  for  the 
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Figure  ^4-1.  Central  Brillouin  zone  for  the  square  lattice. 
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Figure  4-2.  Seven  representative  points  in  the  irreducible  triangle  of 
the  central  Brillouin  zone  for  the  square  lattice. 
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Figure  ^-3.  Central  Brillouin  zone  for  the  hexaconal  lattice. 
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Figure  h-k.    Four  representative  points  of  the  irreducible  triangle  of 
the  central  Briliouin  zone  for  the  hexagonal  lattice. 
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cohesive  energies  per  particle  versus  the  lattice  constant  are  depicted 
in  Figure  ^S.       Initially  these  results  were  interpreted  to  imply  that 
the  square  lattice  structure  possessed  a  lower  equilibrium  cohesive 
energy  per  particle  (-0.035  Hartrees  versus  -0.005  Hartrees)  and  thus 
a  more  stable  equilibrium  geometry  relative  to  the  hexagonal  structure. 

The  effects  of  the  choice  of  a  were  tested  by  first  using  the  molec- 
ular LCAO-Xa  program  to  find  the  value  of  a  which  yields  a  total  energy  of 
atomic  hydrogen  equal  to  -0.5  Hartrees  using  Basis  1.   Further  calculations 
were  performed  using  this  value  of  a  (equal  to  O.78I508)  and  it  was  found 
that  the  cohesive  energy  curve  underwent  only  minimal  changes,  with  the 
equilibrium  cohesive  energy  per  particle  of  the  square  structure  changing 
by  less  than  0.001  Hartrees  with  respect  to  the  equilibrium  cohesive  energy 
per  particle  of  0.035  Hartrees  for  the  square  lattice  calculations  with 
a  equal  to  2/3. 

More  extensive  calculations  were  performed  using  a  (3s, Ip)  basis  set 
for  the  orbital  basis  set  where  p-type  functions  with  exponents  equal  to 
1.0  were  added  to  Basis  1  as  polarization  functions.   This  basis  set 
(henceforth  referred  to  as  Basis  2)  is  presented  in  Table  k-]    along  with 
Basis  1.   In  addition  to  the  slight  improvement  in  the  basis  set,  a  finer 
grid  of  wave  vectors  was  used  in  the  Brillouin  zone,  increasing  the  number 
of  non-equivalent  (by  translation  symmetry)  points  in  the  central  Bril- 
louin to  256.   This  increases  the  number  of  non-equivalent  points,  when 
point  group  symmetry  is  considered,  in  the  Brillouin  zone  from  k   points 
to  30  points  for  the  hexagonal  lattice  and  from  7  points  to  kS   points  for 
the  square  lattice.   Figures  k-G   and  k-J   display  these  two  sets  of  points 
in  their  irreducible  regions  of  the  Brillouin  zone  for  the  square  and 
hexagonal  lattices  respectively. 
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square  (Q)  and  hexagonal  (A)  lattices. 
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Figure  A-6.  Forty-five  representative  points  of  the  irreducible  triangle 
of  the  Brillouin  zone  for  the  square  lattice. 
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Figure  4-7.  Thirty  representative  points  of  the  irreducible  triangle  of 
the  Srillouin  zone  for  the  hexaoonal  lattice. 
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The  resulting  cohesive  energies  per  particle  and  virial  ratios  are 
presented  as  a  function  of  the  lattice  constant  for  the  square  and  hex- 
agonal lattices  in  Tables  1^-2   and  h-3    respectively.   A  graphical  display 
of  the  cohesive  energies  per  particle  is  made  in  Figure  h-8.      The  cohe- 
sive energies  are  calculated  as  the  total  energy  minus  the  energy  of  the 
atomic  monolayer  at  a  lattice  separation  of  100  a.u.   As  can  be  easily 
seen,  these  resul ts  differ  markedly  from  the  results  of  our  preliminary 
calculations.   Our  earlier  inference  that  the  square  lattice  is  the  ener- 
getically preferred  geometry  at  equilibrium  is  obviously  refuted  by  these 
more  recent  results.   In  fact  by  plotting  the  cohesive  energy  per  particle 
versus  the  square  root  of  the  area  of  the  cross  section  in  the  periodic 
plane  of  the  parallelepiped  unit  cell  defined  in  Section  2-1  (thus  giving 
the  cohesive  energy  per  particle  as  a  function  of  the  density  of  the 
lattice  sites)  as  displayed  in  Figure  '♦-9,  we  see  that  the  results  for 
the  two  different  geometries  are  almost  indistinguishable.   Most  of  the 
difference  between  the  two  sets  of  results  may  be  attributed  to  the  in- 
crease in  the  number  of  Brillouin  zone  points,  since  the  improvement  in 
the  basis  from  a  (3s)  basis  to  a  (3s, Ip)  basis  is  only  a  modest  improve- 
ment compared  to  the  increase  of  the  total  number  of  points  in  the  Bril- 
louin zone  from  16  to  256. 

It  is  to  be  noted  that  both  Table  h-2   and  Table  ^4-3  contain  two  sets 
of  results  for  the  lattice  constant  of  5-0  a.u.   Using  the  spin-polarized 
version  of  the  computer  code  at  this  particular  lattice  separation  (and 
only  at  this  particular  choice)  leads  to  two  stable  (in  the  sel f -cons i stent 
iterative  procedure)  solutions  that  correspond  to  a  non-spin-polarized 
(NSP)  and  completely  sp i n-polar i zed  (SP)  solutions.   V/hich  result  is  ob- 
tained depends  on  the  initial  choice  of  fitting  coefficients  to  generate 
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Table  A-2.  Calculated  results  for  the  square  atomic  hydrogen  monolayer. 


Latt  ice 

Separat  ion 

Total  Energy 

Binding 

Energy 

(a.u.) 

(Hartrees) 

-V/T 

(Hartrees) 

(eV) 

2.25 

-0.^883 

1.867 

-0.0354 

-0.963 

2.50 

-0.i*921 

1.995 

-0.0391 

-1.065 

2,65 

-0.4918 

2.062 

-0.0389 

-1.057 

2.75 

-0.4909 

2.103 

-0.0380 

-1.033 

3.00 

-0.4868 

2.188 

-0.0338 

-0.920 

5.003 

-0.4327 

2.187 

+0.0202 

+0.550 

s.oo^^ 

-0.4521 

1.934 

+0.0009 

+0.024 

100.00 

-0.4530 

1.972 

non-spin-polarized  solution 
spin-polarized  solution 
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Table  h-3.    Calculated  results  for  the  hexagonal  atomic  hydrogen  monolayer, 


Lattice 

Separation 

Total  Energy 

Binding 

Energy 

(a.u.) 

(Hartrees) 

-V/T 

(Hartrees) 

(eV) 

2.25 

-0.^823 

1.784 

-0.0293 

-0.797 

2.50 

-0.4917 

1.919 

-0.0332 

-1.040 

2.65 

-0.i»92'i 

1.990 

-0.0394 

-1.073 

2.75 

-0.4922 

2.033 

-0.0392 

-1.066 

3.00 

-0.4892 

2.131 

-0.0363 

-0.987 

5.00^ 

-0.4378 

2.229 

+0.0151 

+0.412 

5.00'^ 

-0.4503 

1.923 

+0.0027 

+0.072 

100.00 

-0.4530 

1.971 

non-spin-polarized  solution 
spin-polarized  solution 
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Figure  k-Q.    Cohesive  energy  per  particle  versus  lattice  constant  for 
square  (□)  and  hexagonal  (A)  lattices. 
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Figure  ^-3.  Cohesive  energy  per  particle  versus  nodified  lattice  constant 
for  square  (Q)  and  hexagonal  (A)  lattices. 
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the  effective  potential,  since  using  fitting  coefficients  corresponding 
to  an  SP  system  of  overlapping  atomic  potentials  yields  the  SP  result, 
while  starting  with  coefficients  from  the  NSP  state  at  3.0  a.u.  yields 
the  NSP  result.   For  the  NSP  results  in  general,  it  was  necessary  to 
rerun  the  calculation  in  the  NSP  mode  of  the  computer  code  in  order  to 
confirm  this  assertion,  since  technical  details  involving  choosing  occu- 
pation numbers  at  wave  vectors  near  the  Fermi  surface  prevented  the  SP 
version  from  converging  on  a  totally  NSP  result.   The  NSP  results  cal- 
culated in  both  the  NSP  and  SP  modes  of  the  computer  code  differed  by 
less  than  10"5  Hartrees  in  the  total  energy.   The  NSP  and  SP  results  at 
a  lattice  constant  equal  to  5-0  a.u.  evidently  indicate  that  the  bands 
corresponding  to  the  different  spins  collapse  at  a  lattice  constant 
slightly  less  than  5.0  a.u.   This  contention  was  tested  by  calculating 
the  cohesive  energies  per  particle  for  both  the  NSP  and  SP  modes  at 
lattice  constants  equal  to  6.0  a.u.,  8.0  a.u.,  and  12.0  a.u.   Figure 
k-]0    illustrates  the  results  for  these  calculations  and  demonstrates  the 
apparent  crossing  of  the  NSP  state  and  the  SP  state  at  a  lattice  constant 
between  4.0  a.u.  and  5.0  a.u.   A  similar  phencmenon  has  been  reported 
for  the  crystalline  vanadium  system  (Hattox,  Conklin,  Slater,  and  Trickey, 
1973)  which  exhibits  the  same  collapse  of  bands  of  different  spins  to  a 
doubly  degenerate  band  as  the  lattice  constant  is  decreased  from  the 
atomic  1 imi  t. 
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Figure  ^-10.  Spin-polarized  and  non-spin-polarized  solution  cohesive 

energies  per  particle  versus  lattice  constant  for  the  square 
lattice.   Dashed  portion  of  spin-polarized  curve  is 
extrapolated  from  the  solid  portion  of  the  line. 


65 


^-2.  Atomic  Beryllium  Monolayer 

Beryllium  is  an  interesting  system  to  study  using  the  metiiods  des- 
cribed in  this  work  for  several  reasons.   Recent  experimental  measurements 
of  the  work  function  of  the  (OOOl)  surface  of  beryllium  by  Green  and 
Bauer  (1978)  indicate  that  beryllium  has  a  work  function  of  5-1  eV ,  a 
value  much  higher  than  previous  measurements  had  indicated.   Until  the 
results  of  Green  and  Bauer  the  recommended  value  of  the  work  function 
according  to  the  standard  reference  (Fomenko,  1966)  was  only  3-92  eV. 
Although  Green  and  Bauer  attribute  this  discrepancy  in  the  work  function 
to  possible  oxidation  of  the  beryllium  surface  in  the  earlier  work,  these 
differences  in  experimental  results  pose  a  question  worthy  of  theoretical 
investigation.   The  Interpretation  of  certain  peaks  in  the  Auger  spectra 
of  beryllium  (Suleman  and  Pattlnson,  1971)  also  poses  questions  appropriate 
for  theoretical  investigation  using  surface  electronic  structure  computa- 
tional methods. 

In  addition  to  the  questions  resulting  from  experimental  study,  beryl- 
lium Is  an  interesting  prototype  system  for  studying  the  surface  of  metals. 
With  only  four  electrons  per  atom,  beryllium  should  require  relatively 
modest  amounts  of  computational  effort  for  an  adequate  treatment  of  the 
surface.   Also  beryllium  is  only  v;eak1y  bound  as  a  dimer,  with  a  dissoci- 
ation energy  which  is  experimentally  estimated  to  be  about  0.7  eV  (Gaydon, 
1968).   The  weak  binding  of  the  dimer  is  also  indicated  by  theoretical 
calculations  such  as  the  configuration  interaction  calculations  by  Bender 
and  Davidson  (1967),  which  yield  a  strictly  repulsive  ground  state  poten- 
tial energy  curve,  and  the  sel f -cons  I  stent  electron  pair/coupled  electron 
pair  approximation  (SCEP/CEPA)  calculations  by  Dykstra,  Schaefer,  and 
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Meyer  (1977)  which  leads  them  to  estimate  the  calculated  dissociation 
energy  to  be  about  0.03  eV  at  an  equilibrium  separation  of  about  8.5  bohr. 
Hartree-Fock  calculations  by  Bauschl icher ,  Liskow,  Bender,  and  Schaefer 
(1975)  and  by  Jordan  and  Simons  (1977)  indicate  the  smallest  cluster  that 
exhibits  appreciable  binding  (0.'46  eV  per  atom)  contains  four  beryllium 
atoms.   This  indicates  that  the  cohesive  energy  of  the  solid  (and  pre- 
sumably the  monolayer)  may  be  primarily  due  to  the  de loca 1 i za t ion  of  the 
electron  charge  over  a  collection  of  sites.   The  full  effect  of  this  de- 
localization  will  be  investigated  by  comparing  cohesive  energies  of  the 
infinitely  extended  monolayer  using  the  methods  described  in  this  work 
with  the  binding  energies  calculated  by  Bauschlicher  (1976)  on  beryllium 
clusters  and  with  experimental  results  for  the  real  surface. 

We  have  performed  calculations  on  the  atomic  beryllium  monolayer 
using  the  (6s,  2p)  basis  given  in  Table  '4-'4.   The  6  s-type  functions  are 
taken  from  van  Diujneveldt ' s  (1971)  6s  basis  for  beryllium.   The  2  p-type 
functions  are  the  same  as  used  by  Bauschlicher  (1976)  in  his  cluster  cal- 
culations where  the  p-type  functions  were  used  with  van  Du i jneveldt ' s  9s 
basis  for  beryllium.   The  p-type  functions  are  introduced  as  polarization 
functions,  since  no  p  orbitals  will  be  occupied  in  the  ground  state  of  the 
beryllium  atom.   However  since  beryllium  contains  the  greatest  number  of 
electrons  of  the  elements  that  have  only  s  orbitals  occupied,  we  may  antic- 
ipate that  the  p  functions  will  have  a  greater  contribution  than  equivalent 
calculations  for  hydrogen,  helium,  or  lithium.   Thus  questions  of  the  bal- 
ance between  the  s  functions  and  p  functions  should  be  considered.   In 
his  study  of  basis  set  effects,  van  Diujneveldt  (1971)  states  that  s/p 
ratios  greater  than  5/2  (for  2  p-type  functions)  lead  to  unbalanced 
behavior  in  Hartree-Fock  calculations  on  the  nitrogen  dimer.   Since 
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Table  h-k .    Basis  sets  for  the  beryllium  monolayer. 


Orbital  basis 


Charge  fitting 

bas  is 


Exchange  fitting 
bas  i  s 


Type 


s 
s 
s 
s 
s 
s 

X 
X 

y 
y 

z 
z 


Exponent 

0.067376 
0.198210 
1.767558 
6.819528 
30.827597 
204.90614'* 
0.509 
0.118 
0.509 
0.118 
0.509 
0.118 


Type 


s 
s 
s 
s 
s 
s 
s 
s 

x2+y2 
x2+y2 

z2 


Exponent 

0.134752 

0.396420 

3.535116 
13.639056 
61.655194 
409.81229 

1.018 

0.236 

1.018 

0.236 

1.018 

0.236 


Type 


s 
s 
s 
s 
s 
s 
s 


x2+y2 
x2+y2 
22 


Exponent 

0.0449173 
0.132140 
1.178372 
4.546352 
20.551731 
136.604096 
0.33933 
0.079667 
0.33933 
0.079667 
0.33933 
0.079667 
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beryllium  has  no  occupied  p  orbitals  in  the  atomic  ground  state  (while 
nitrogen  does)  we  feel  that  our  6s/2p  ratio  is  acceptable,  especially  in 
the  light  of  the  results  by  Bauschlicher  (1976)  that  a  3s/lp   ratio  yields 
reasonable  results.   Preliminary  calculations  on  the  dimer  were  performed 
with  the  above  basis  using  the  molecular  LCAO-Xa  program  (Mintmire,  1979) 
and  a  plot  of  the  binding  energy  versus  the  nuclear  separation  is  presented 
in  Figure  '4-II.   These  results  are  consistent  with  the  previously  men- 
tioned results  of  Dykstra,  Schaefer,  and  Meyer  and  demonstrate  the  weak 
binding  of  the  beryllium  dimer  as  well  as  provide  a  check  on  the  balance 
of  the  (6s,  2p)  basis  set. 

Since  beryllium  possesses  a  hexagonal  close  packed  structure  in 
the  solid,  the  calculations  were  performed  on  a  hexagonal  lattice  of  beryl- 
lium atoms.   This  structure  is  equivalent  to  a  single  layer  from  the  (0001) 
surface  of  the  crystalline  solid.   The  beryllium  lattice  possesses  a  Bril- 
louin  zone  with  the  same  symmetry  as  that  of  the  hexagonal  atonic  hydrogen 
monolayer  previously  displayed  in  Figure  k-},.      We  used  the  same  array  of 
30  non-equivalent  points  in  the  Brillouin  zone  for  the  beryllium  as  for 
the  hexagonal  atomic  hydrogen  monolayer  calculations. 

Table  h-S   presents  some  results  of  our  calculations  giving  the  co- 
hesive energies  and  virial  ratios  of  the  beryllium  monolayer  as  a  function 
of  the  lattice  spacing.   These  results  and  Figure  ^4-12  imply  that  the 
beryllium  monolayer  has  an  equilibrium  lattice  separation  of  '♦.I  bohrs  and 
an  equilibrium  cohesive  energy  of  -2.6^4  eV  per  particle.   Since  the  virial 
theorem  is  valid  for  the  Xa  method  (Slater,  1972b),  the  virial  ratio  -V/T 
should  equal  2  at  the  minimum  of  the  cohesive  energy  curve  where  the  de- 
rivative of  the  cohesive  energy  with  respect  to  the  lattice  spacing  equals 
zero.   We  see  in  Table  h-S    that  the  location  of  the  cohesive  energy 
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Table  4-5. 

Calculated  results 

for  the 

atomic  beryl  1  iuir 

1  monolayer. 

Lattice 

Separation 

Total  Energy 

Bindi 

ng  Energy 

(a.u.) 

(Hartrees) 
-14.293628 

-V/T 

(Hartrees) 
-0.092185 

(eV) 

3.80 

1.99370 

-2.5084 

3.90 

-14.296508 

1.99553 

-0.095065 

-2.5867 

4.00 

-14.298019 

1.99880 

-0.096576 

-2.6278 

4.10 

-14.298351 

2.00149 

-0.096908 

-2.6369 

4.25 

-14.297454 

2.00586 

-0.096011 

-2.6125 

4.321 

-14.296885 

2.00566 

-0.095442 

-2.5970 

4.50 

-14.293275 

2.00837 

-0.091832 

-2.4987 

5.00 

-14.275816 

2.01493 

-0.074373 

-2.0237 

6.00 

-14.234844 

2.01263 

-0.033401 

-0.9088 

100.00 

-14.201443 

1.99901 
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minimum  (predicted  by  interpolation  of  the  virial  ratios  to  yield  a  pre- 
dicted lattice  constant  of  ^.05  a.u.  and  an  estimated  cohesive  energy  per 
particle  at  that  lattice  constant  of  -2.63  eV)  is  in  close  agreement  with 
that  predicted  by  the  actual  cohesive  energy  curve.   These  equilibrium 
results  differ  slightly  from  the  bulk  equilibrium  separation  of  ^4.321 
bohrs  and  equilibrium  cohesive  energy  of  -3-32  eV  per  particle  (Kittel, 
1976),  but  not  more  than  could  be  attributed  to  the  difference  between 
monolayer  and  bulk  solid  beryllium.   Bauschlicher  (1976)  has  performed 
calculations  using  the  Hartree-Fock  method  on  various  size  beryllium  clus- 
ters with  the  same  internuclear  separation  as  the  bulk  solid.   For  the 
various  monolayer  clusters,  Bauschlicher  reports  binding  energies  of  +0.22 
eV  per  particle  for  a  three  atom  cluster,  +0.07  eV  per  particle  for  a  six 
atom  cluster,  -O.I6  eV  per  particle  for  a  seven  atom  cluster,  and  -0.63  eV 
per  particle  for  a  fourteen  atom  cluster.   The  largest  cluster  Bauschlicher 
considered  was  a  twenty-two  atom  cluster  in  tv^o  layers  which  yielded  a 
computed  binding  energy  of  -O.89  eV  per  particle.   These  cluster  results 
apparently  indicate  a  trend  in  binding  energies  as  the  cluster  size  in- 
creases which  should  be  consistent  with  our  equilibrium  cohesive  energy 
of  -2.6^  eV  per  particle,  although  one  may  infer  that  much  larger  clusters 
are  necessary  for  a  proper  treatment  of  the  cohesive  energy  of  the  beryl- 
lium monolayer.   This  Is  a  point  in  favor  of  our  earlier  contention  that 
a  primary  cause  of  the  binding  energy  of  the  beryllium  monolayer  Is  the 
del ocal I zat Ion  of  the  electronic  charge  over  many  nuclear  sites. 

Band  energies  were  Interpolated  for  ^4096  points  in  the  Brlllouln  zone 
using  a  scheme  adapted  from  the  method  proposed  by  Monkhorst  and  Pack 
(1976)  for  three-dimensional  Brillouin  zone  interpolations.   This  method 
uses  a  discrete  Fourier  transform  over  the  256  evenly  spaced  points  In  the 
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complete  Brillouin  zone  for  which  the  band  energies  are   computed.   Band 
energies  at  intermediate  points  were  then  evaluated  using  the  resulting 
Fourier  expansion.   Figure  A-13  displays  the  valence  bands  and  lower  empty 
bands  computed  for  the  beryllium  monolayer  at  a  lattice  spacing  of  ^.10  a.u. 
The  shape  of  the  bands  resemble  qualitatively  the  bands  reported  by  Loucks 
and  Cutler  (196^*)  for  crystalline  beryllium,  if  we  compare  our  results 
with  the  solid  energy  bands  for  wave  vectors  lying  in  planes  parallel  to 
the  plane  of  periodicity  of  the  monolayer.   One  may  assume  that  the  dis- 
crepancy between  our  bands  and  those  of  Loucks  and  Cutler  is  due  to  the 
differences  between  the  monolayer  and  the  crystal,  and  that  the  addition 
of  more  layers  would  lead  to  multilayer  bands  more  in  agreement  with  those 
of  the  crystalline  beryllium  solid.   Further  work  is  currently  progressing 
in  the  direction  of  calculations  involving  more  than  one  layer  of  beryllium. 

We  may  estimate  the  work  function  of  beryllium  from  the  negative  of 
the  Fermi  energy  of  our  calculations  (Kittel,  1976)-   V/e  find  that  our 
Fermi  energy  for  the  monolayer  of  -3.80  eV  agrees  quite  well  with  the 
recommended  value  of  Fomenko  (I966).   However  the  discussion  by  Trickey 
and  Worth  (1977)  indicates  that  while  the  choice  of  a   equal  to  2/3  will 
yield  reasonable  cohesive  energies  and  band  widths,  one  expects  that  this 
value  of  a  will  lead  to  underestimating  band  gaps  and  the  Fermi  energy. 
Thus  the  fact  that  we  anticipate  our  Fermi  energy  to  underestimate  the 
work  function,  and  that  increasing  the  number  of  layers  will  probably 
change  the  Fermi  energy  in  a  manner  not  easily  predictable,  our  estimated 
work  function  of  3.8O  eV  is  not  necessarily  in  conflict  with  the  reported 
experimental  work  function  of  5.10  eV  (Green  and  Bauer,  1978). 

The  density  of  states  for  the  beryllium  monolayer  was  computed  using 
the  interpolated  set  of  band  energies  and  is  presented  in  Figure  'l-H. 
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Figure  ^-13.  Band  energies  for  the  hexagonal  beryllium  monolayer.   The 
core  band  located  at  -100.93  eV  is  omitted. 
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Figure  k-]k.    Total  density  of  states  for  the  hexagonal  beryllium  mono- 
layer.  The  solid  line  shown  in  the  inset  is  the  total 
density  of  states  computed  for  crystalline  beryllium  by 
Loucks  and  Cutler  (I964) . 
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The  results  are  also  consistent  and  in  qualitative  agreement  with  the 
density  of  states  computed  for  the  crystalline  beryllium  solid  (Loucks 
and  Cutler,  \S6k) .      The  relatively  flat  region  in  the  density  of  states 
at  energies  around  -10.  eV  in  Figure  h-\k    is  understandable  when  one  no- 
tices that  the  bands  in  Figure  ^-13  are  nearly  free  electron  like  in  char- 
acter.  One  can  demonstrate  that  the  density  of  states  for  a  two-dimen- 
sional electron  gas  will  be  a  step  function.   Thus  the  flat  region  of  our 
density  of  states  is  due  to  the  nearly  free  electron  like  behavior  of  the 
bottom  of  the  lowest  band  in  Figure  '4-12. 

Musket  and  Fortner  (1971)  have  reported  Auger  Electron  Spectra  (AES) 
results  for  beryllium  Indicating  two  peaks  at  92  eV  and  10^4  eV.   Later 
work  by  Suleman  and  Pattinson  (1971)  indicates  that  the  AES  peak  at  92  eV 
was  evidently  due  to  impurities  on  the  beryllium  surface,  although  they 
agreed  with  the  earlier  conclusion  of  Musket  and  Fortner  that  the  peak  at 
104  eV  was  due  to  a  (2p,  2p)  KVV  transition  at  the  solid  surface.   The 
energy  of  this  transition  will  equal  the  energy  difference  of  an  electron 
in  the  2p  band  (e   )  and  the  energy  of  the  hole  it  fills  in  the  core  K 
shell  (e^)  minus  the  energy  required  to  excite  an  electron  from  some  dif- 
ferent  location  in  the  2p  band  {c'    ).   We  may  thus  find  an  estimate  for 
the  peak  of  the  above  mentioned  transition  by  taking  the  difference 

Ae  =  e„  +  e-i   "  El/-   To  find  appropriate  bounds  for  our  estimate  (since 
Zp     /p     K 

we  do  not  anticipate  monolayer  results  more  accurate  than  2  or  3  eV  in  the 
light  of  our  previous  discussion  of  the  effects  of  adding  more  layers  to 
the  beryllium  system)  we  may  consider  e-   and  e'      jointly  equal  either  to 
the  bottom  of  the  second  band  in  Figure  4-13,  which  corresponds  to  the 
2p   band,  or  to  the  Fermi  energy.   Since  we  anticipate  that  the  band  energy 
for  the  core  band  will  not  be  very  accurate,  we  will  use  the  reported 
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results  of  -111  eV  plus  the  work  function  (Siegbahn,  et  al.,  196?). 
Depending  upon  whether  we  use  the  experimental  work  function  or  our  esti- 
mated work  function,  we  find  a  range  of  105  to  109  eV  for  our  estimate 
of  the  (2p,  2p)  KVV  Auger  peak.   Thus  we  find  one  more  source  of  agree- 
ment of  our  monolayer  results  with  those  of  the  bulk  solid,  indicating 
that  even  these  simple  monolayer  calculations  provide  useful  information 
for  describing  the  surface  of  the  bulk  solid. 

These  calculations  were  performed  on  an  Amdahl  V?  computer  and  re- 
quired approximately  20  minutes  of  CPU  time  for  each  value  of  the  lattice 
constant.   Preliminary  calculations  have  been  attempted  for  the  beryllium 
double  layer  with  relaxed  precision  tolerances  in  order  to  reduce  the  CPU 
time  required.   Unfortunately  these  preliminary  calculations  required 
approximately  k   hours  of  CPU  time  per  lattice  constant  for  the  double 
layer  and  yielded  unsatisfactory  results  due  to  precision  problems  involv- 
ing the  more  diffuse  basis  functions.   It  is  felt  that  the  inclusion  of 
routines  to  incorporate  the  mixed  lattice  expansion  algorithm  for  evalu- 
ating the  V.^.|^(k.^^)  and  X  .|^(k_^^)  integrals  described  in  Chapter  ill  will 
both  reduce  the  total  required  computational  time  and  increase  the  pre- 
cision of  the  evaluated  integrals  involving  diffuse  functions.   Steps  are 
currently  being  taken  to  incorporate  these  techniques  into  the  computer 
code  in  order  to  perform  calculations  on  the  beryllium  double  layer  and 
other  more  complex  two-d Imens ional 1 y  periodic  systems. 


CHAPTER  V 
SUMMARY 


The  primary  purpose  of  this  work  has  been  to  develop  and  test  a  new 
approach  for  calculating  the  electronic  structure  of  thin  films,  with  the 
ultimate  objective  of  studying  the  electronic  properties  of  solid  sur- 
faces.  The  calculations  presented  in  this  dissertation  hopefully  have 
indicated  some  of  the  possible  applications  of  our  approach,  although  a 
complete  understanding  of  the  strengths  and  ^^/eaknesses  of  the  approach 
outlined  in  this  dissertation  will  require  much  more  extensive  testing. 
The  basic  approach  discussed  herein  has  relied  on  the  following  assump- 
tions: 

1)  Properties  characteristic  of  the  surface  can  be  reliably 
studied  using  a  thin  layer  model; 

2)  the  coupling  of  nuclear  motion  and  electronic  motion  may  be 
neglected  so  that  electronic  properties  can  be  calculated 
reliably  within  the  Born-Oppenhe imer  approximation;  and 

3)  the  local  density  functional  (Xa)  formalism  yields  results 
which  are   useful  in  understanding  the  electronic  properties 
of  surfaces. 

The  first  assumption  is  intuitively  reasonable  given  sufficiently 
many  layers  in  the  film,  and  its  use  is  supported  by  work  on  d-band 
metals  which  indicates  that  the  surface  effects  exhibited  by  films  15  to 
to  20  layers  thick  accurately  simulate  the  effects  of  the  bulk  surface 
(Cooper,  1977).   The  Born-Oppenheimer  approximation  is  a  common  one  in 
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both  molecular  and  solid  state  calculations,  and  barring  evidence  to  the 
contrary  for  a  particular  case,  we  may  assume  this  approximation  to  be 
reasonable.   The  third  assumption  has  been  shown  empirically  to  yield 
results  in  reasonable  agreement  with  experiment  in  a  multitude  of  cases 
(for  a  review  of  applications  of  the  Xa  formalism  see  Slater,  197!,,  and 
Connolly,  1976).   The  final  test  of  all  three  assumptions,  of  course,  is 
that  of  comparison  of  results  for  a  series  of  different  systems  using 
the  above  model  with  the  experimental  results  for  those  systems. 

In  addition  to  the  previously  mentioned  assumptions  used  in  estab- 
Mshing  the  mathematical  model,  various  approximations  have  been  made  in 
finding  a  computational  solution  for  the  model.   These  approximations, 
which  are  capable  of  refinement  to  any  limit  of  precision  (or  at  least  to 
the  limit  of  precision  of  the  computational  equipment)  with  a  correspond- 
mg  increase  in  computational  time,  fall  in  four  basic  categories: 

1)  the  representation  of  orbitals  with  a  finite  basis  (the  LCAO 
approximation) ; 

2)  the  general  computational  approximation  of  evaluating  quantities 
to  within  some  predetermined  tolerance; 

3)  the  procedure  for  fitting  the  charge  density  and  the  cube  root 
of  the  charge  density;  and 

A)   the  approximate  treatment  of  the  Brillouin  zone  introduced  through 
the  use  of  periodic  boundary  conditions. 
Of  these  approximations,  the  first  two  are  well  known  and  no  detailed  dis- 
cussion of  these  points  will  be  made.   The  effects  of  finite  basis  sets 
are  discussed  by  Schaefer  (1972,  pp.  56-83)  while  the  effects  of  precision 
tolerances  are  discussed  by  dementi  and  Mehl  (1971).   No  attempts  have 
been  made  to  optimize  basis  sets,  although  bases  optimized  for  atomic 


80 


hydrogen  and  beryllium  in  other  works  have  been  used  in  our  reported 
calculations.   Appropriate  tolerances  were  found  so  that  the  total  ener- 
gies were  calculated  to  an  estimated  precision  of  10   Hartrees.   The 
effects  of  using  different  numbers  of  points  in  the  Brillouin  zone  has 
been  examined  briefly  for  the  hydrogen  monolayer,  with  the  conclusion 
that  an  insufficient  number  of  points  can  yield  radically  different  co- 
hesive energies  from  results  using  an  adequate  number  of  Brillouin  zone 
points.   Examining  the  difference  of  band  energies  between  adjacent  points 
in  the  Brillouin  zone  for  the  beryllium  calculations  using  256  points  in 
the  Brillouin  zone  yields  an  upper  bound  to  the  error  in  the  total  energy 
of  about  0.05  Hartrees  for  this  particular  source  of  error.   No  testing 
has  been  made  of  the  errors  introduced  by  the  charge  density  fitting 
procedures,  although  we  anticipate  that  errors  in  the  total  energy  due  to 
this  source  of  error  will  be  less  than  0.001  Hartrees  if  the  fitting 
scheme  for  thin  films  behaves  similarly  to  the  fitting  scheme  for  the 
molecular  case  (Dunlap,  Connolly,  and  Sabin,  1979a). 

The  physically  reasonable  results  produced  for  the  hydrogen  mono- 
layer indicate  that  there  are  no  inherent  flaws  in  the  mathematical  for- 
malism.  The  beryllium  monolayer  results  further  indicate  the  feasibility 
of  our  approach  as  well  as  demonstrating  the  applicability  to  experimen- 
tally interesting  systems.   The  comparison  of  our  theoretically  estimated 
values  for  the  equilibrium  lattice  spacing,  equilibrium  cohesive  energy 
per  particle,  and  work  function  with  the  experimentally  determined  quan- 
tities for  the  beryllium  system  and  with  other  calculated  results  indicate 
that  even  such  a  simple  model  as  the  monolayer  yields  useful  results. 
The  primary  difficulty  with  our  approach  at  this  time  is  the  amount  of 
computer  time  required  to  produce  adequately  precise  results  for  systems 
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larger  than  the  atomic  beryllium  monolayer.   Current  work  to  incorporate 
the  reciprocal  lattice  expansion  techniques  described  in  Appendix  C  as 
well  as  structural  changes  in  the  computer  code  to  optimize  the  order  of 
operation  should  lead  to  a  substantial  reduction  in  required  computer 
time. 

In  summary  we  may  state  that  this  dissertation  hopefully  represents 
the  beginning  of  an  extended  study  of  the  electronic  properties  of  sur- 
faces.  Needless  to  say,  the  methods  presented  herein  are  not  a  finished 
product,  but  rather  the  first  stage  of  development  of  what  we  believe  to 
be  a  line  of  study  that  can  be  extremely  fruitful  in  terms  of  understanding 
the  basic  properties  of  solid  surfaces. 


APPENDIX  A 
HERMITE-GAUSSIAN  FUNCTIONS  AND  INTEGRALS 


A-1.  Definition  of  Hermi te-Gauss ian  Functions 

The  basis  functions  used  in  this  work  are  functions  which  are  com- 
posed of  specific  combinations  of  Cartesian  Gaussians  (which  are  of  the 
form  X  y  z  exp(-ar  },  where  i,  m,  and  n  are  non-negative  integers). 
These  functions  may  be  referred  to  as  Hermi te-Gauss ian  functions,  which 
is  the  nomenclature  introduced  by  Zivkovic  and  Maksic  (1968)  in  their 
discussion  of  the  use  of  these  functions  for  molecular  calculations.  A 
Hermi  te-Gaussian  function  f  (n ,  a,  A_;  r)     is  defined: 

f(n,  a,  A;  0   =   3^  exp(-a  |_r  -  Ap)  (A-l) 


here  n  refers  to  a  set  of  three  non-negative  integers  (n^,  n2,  nj)  which 


w 


define  the  operation  8.: 


n             a"l    a"2    ^"3 
a"  =   3 9 3 (A-2) 

^     8a"1   3A"2   3a"3 


This  notation  is  consistent  with  that  of  Zivkovic  and  Maksic,  although 
other  workers  (Golebiev'^ski  and  Mrozek,  1973;  McMurchie  and  Davidson, 
1978)  have  introduced  slightly  different  notations.   The  term  "Hermite- 
Gaussian  functions"  derives  from  the  fact  that  given  the  form  of  the 
Hermite  polynomials: 
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H  (x)   =   (-1)"  exp(x^)  a"  exp(-x^)  (A-3) 

n  A. 

the  Hermi  te-Gauss  ian  function  f(n,  a,  A_;  r)    may  be  expressed: 

f(n,  a,  A;  r)  =  a^^^  H   [a^/^(x-A  )]  H   [a^^^y-A  )] 
—  —         n  1        X    n2        y 

H   [a'/^(2-A  )]  exp(-alr  -  A|^)     (A-^) 

where  n  denotes  the  sum  of  the  components  of  n  (i.e.,  n  =  nj  +  n2  +  n3). 
These  functions  are  useful  in  computation  since  they  enable  us  to  per- 
form the  calculation  of  integrals  and  other  numerical  data,  as  pointed 
out  by  Zivkovic  and  Maksic  (1968),  in  three  steps: 

1)  for  any  integration  or  summation  which  does  not  involve  A^  and 
is  linear  in  f(n,  a,  A^;  r_)  ,  bring  the  differentiation  operator 
outside  the  integration  or  summation; 

2)  calculate  the  appropriate  integral  or  sum  for  spherically  sym- 
metric Gaussians  in  closed  form;  and 

3)  apply  the  differentiation  operator  to  the  analytic  form  of  the 
resultant  expression. 

As  an  example  we  may  compute  the  overlap  integral  between  two 
Hermi te-Gauss ian  functions  centered  on  two  different  centers  A  and  B. 

d^r  f(n,a,A;  r)    f(n',b,B^;  _r) 

=   9^  3^'  jd^r   exp[-a|r-A|^  -  b|£-B_|2]         (A-5) 

=  a^  3^'  [Ti/(a+b)]^^^  exp[-  ab|A  -  B|^/(a+b)] 
-  (-1)"'  [TT/(a+b)]^^^-  f(n-i-n',  ab/(a+b),  A-B;  0) 
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A-2.  Products  of  Two  Hermi te-Gauss ian  Functions 

The  overlap  integrals  of  Equation  A-5  are  a  simple  example  of  inte- 
grals which  involve  the  product  of  two  Hermi te-Gauss ian  functions.   Such 
integrals  may  be  evaluated  with  less  complexity  by  decomposing  the  prod- 
uct of  two  Hermi te-Gaussian  functions  into  two  factors.   Consider  the 
product  of  two  Hermi te-Gaussian  functions: 


f(n,  a,  A;  r)   f(n',  b,  B;  r) 


=  a^  9^'  exp[-al£  -  Al^]  exp[-blt^  -  BJ^]       (A-6) 


=   3?   9^"    exp[  ^  lA-Bl^]    exp  [- (a+b)  |  r   -  f4_lil:|2] 


Defining  the  quantitities  P^  and  Q: 

P_  =  A  -  B_  (A-7) 

A 

and 

Q  =   (aA  +  bB_)  /  (a  +  b)  (A-8) 

we  may  express  the  product  of  two  Hermi te-Gaussian  functions  in  terms 
of  functions  g  ( P_)  and  h  ( r_  -  Q)  : 

f(n,  a,  A;  r_)  f(n',  b,  B;  r)    =  9^  9^'  g(P)  h  (j^  -  Q)  (A-9) 

where 

g(P)  =  exp[  ^  |P|2]  (A-10) 

and 

h(£  -  Q)  =  exp[  -{a+b)|£  -  Q|^]  (A-ll) 
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Repeated  taking  of  derivatives  yields: 

:   :.  "1    "2   "3   ni'    "2'    "3* 

9^  3    g(P_)  h(£-Q)  =    Z     E     E     E     E     E 

-  -  ki=0  k2=0  k3=0  ki'=0  k2'=0  k3'=0 

A    A       A       A 

X  O^"*"  3g'"'''  h(r-Q)}  (A-12) 


We  may  rearrange  the  sums  over  the  indices  k.  and  k. ' ,  using  the  new 
indices  K.  and  K. ' ,  where 

K.  =  k.  +  k.'  (A-13) 

III  \        J/ 

and 


K.  '  =  k.'  (A-li») 

II  ^  ' 

This  rearrangement  yields  the  expression: 


n      n'  ni+^i'       n2+n2'       nj+nj' 

/3      3g      g(P_)    h(_r-Q_)    =  E  E  E 

-     -  Ki=0  K2=0  K3=0 


min[ni',Ki]  min[n2',K2]  min[n3',K3] 

X  E  E  E 

Ki'=max[0,Ki-ni]      K2 '=max [0,K2-n2]      K3 '=max [O, Ks-ns] 

X      f      "1      1     f      "2      1     f      "3      1     fni'l     (^2']     (^3') 
^Ki-Ki'^     k2-l<2''^     ^K3-K3'-'     ^i'^     ^2 '  ^     ^-Kg'-' 


X      {3^"^'    3^'    g(P)}    {3^"^^^'    3^'"^'    h(r-Q)}  (A-I5) 


Using   the    relationships 


l£    =  is.    =      l£ 

3P.         3A.  3B. 

I  I  I 


(A-16) 
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and 


dh 

a   +  b 

9h 

a   +   b 

9h 

aQ, 

a 

3A. 

1 

b 

9B. 
— 1 

(A-17) 


we  may  rewrite  Equation  A-15: 

^A   ^b'    9^^^    ^^L-1)    =      I  9(K,    n,    n',    a,    b)    {9^  g ( P) } 

K 

X,      r  .n+n  '  -K   ,  /      ^\  , 

Oq  h(£-Q)}  (A- 18) 


AAA 


The   function  g(K,    n,    n',    a,    b)    is   defined: 

'.      -      <s  3  n-K+K'    n'-K'  , 

S(K.    n.   n.,   a.    b)    -     n      (    r      ^  t^,.,      (,%    , )    ("^i,))        (A-19) 

i-i      K..         va+bj  1       1  I 

where  the  indices  K. '  have  as  their  respective  ranges: 

max[0,  K.  -  n.]  <  K. '  <  min[n.',  K.]  (A-20) 


Furthermore  we  may  use  the  relationships: 

A 

3p  g(P.)  =  f(K,  ab/(a+b),  P;  0)  (A-21) 

and 


A       A 

A    A       A 


3q''"'"'^  h(£^-Q)  =  f(n+n'-K,  a+b,  Q;  _r)  (A-22) 

to  yield  the  expression 

f(n,  a,  A;  r_)  f(n',  b,  B^;  r) 

=      Y.     g(K,  n,  n',  a,  b)  f(K,  ab/(a+b),  P^;  0) 
K 

f(n+n'-K,  a+b,  Q;  r)  (A-23) 
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Thus  we  see  that  Equation  A-23  may  be  used  to  reduce  integrals  involving 
the  product  of  two  Hermi te-Gauss ian  functions  of  the  same  coordinate  r 
to  a  sum  of  integrals  involving  a  single  Hermi te-Gauss ian  function  of  the 
coordinate  r_. 

A-3.  Nuclear  Attraction  Integrals 

The  electrostatic  interaction  integral  of  a  charge  distribution 
described  by  a  Hermi  te-Gauss  ian  function  f(n,  a,  A^;  _r_)  with  a  nuclear 
center  of  charge  Z   located  at  the  position  vector  £  may  be  denoted: 


3. 


[f(n,a,A;£)  |  Z^]  =  Z^  /d^r  f(n,a,A;_r_)  i   _  ^i 


(A-2M 


We  may  use  Boys'  (1950)  relation  for  spherically  symmetric  Gaussian  func- 
tions: 

[f(0,a,A;_f)  I  Z^]  =  A  Fo(aT^)  (A-25) 

to  evaluate  this  integral,  where  A  is  the  constant: 


A  =  277  Z^  /  a  (A-26) 


and  F  (s)  is  the  incomplete  gamma  function; 


,1 


Fo(s)  =   /^  du  exp[-su^]  (A-27) 

0 

The  vector  T^  is  defined  as 

T  =  A  -  £  (A-28) 

Thus  the  electrostatic  interaction  integral  in  Equation  A-'ik   may  be 
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expressed: 

[f(n,  a,  A;  r_)  |  Z^]  =  3^  Fo^^"^^^  (A-29) 

McMurchie  and  Davidson  (1978)  have  introduced  the  function  R(n,  a,  r) : 

A 

R(n,  a,  r)    =      a"  F^Car  )  (A-30) 

where  we  have  employed  slightly  different  notation  from  that  of  McMurchie 
and  Davidson.   This  function  may  be  expressed  in  the  expansion: 


o(^      ,      ^^    _      T    ^n-k   „n-2k    /    ,^k      ni-2ki      n2-2k2      n3-2k3    _        /      2v 
R(n,    a,    n    =      E   a  2  (-1)      x^  y  '^   z    ^        ^F   _,  (ar    ) 

/v  n  ~  K 

k 


A  A 


;^l'  "?•'  n^! 

ki!    (ni-2ki)!      k2!    (n2-2k2)  !      ^3'-    (n3-2k3)! 

where  the  incomplete  gamma  functions  F  (s)  may  be  defined 


(A-3I) 


n 


•^     2n    r   2, 


F^(s)  =  /  du  u   exp[-su^]  (A-32) 


and  the  interval  over  which  the  indices  k.  range  is: 

0  <  2k.  <  n.  (A-33) 

I     I 

The  function  R(n,  a,  r)    may  be  evaluated  using  recursion  relationships 

if  we  define  the  quantity  T   ,    .: 

n,l,m,j 

"""n.l.m.j    "    (-/a)"'^'"^"'      /     du      {H^[ux/a]    H^  [uy/a]    H^[uz/a] 
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This  quantity  satisfies  the  following  relationships: 

^O.O.O.j  =  ^-^^^'^  (^-35) 


and 


T  =  R(n,  a,  r)  (A-36) 

ni,n2,n3,0  — 


If  we  use  the  method  described  by  McMurchie  and  Davidson  (1978)  to  eval- 

uate  the  incomplete  gamma  functions  F  (s) ,  the  functions  R(n,  a,  r)  may 
be  evaluated  using  the  recursion  relationships: 

T        ^,  •  =  -2a  {zT         .  ,  +  noT         ,  .^, }  (A-37) 
ni,n2,n3+l,j          ni,n2,n3,j  +  l     "^  n^  ,n2,n3-l  ,  j  +  1 

T     _^,     .  =  -2a  {yT         .  ,  +  noT      ,     .  ,}  (A-38) 

ni,n2+l,n3,j       ^  nj  ,n2  .nj  ,j  +  l     ^  rij  ,n2-l  .nj.j  +  T  ^  ^' 


and 


T  _^,       .  =  -2a  {xT         .,,  +  niT   ,       .  ,}     (A-39) 
ni  +  l,n2,n3,j  ni,n2,n3,j  +  l     ^  ni-1  ,n2,n3  ,j  +  r 


Thus  we  see  that  the  electrostatic  interaction  integral  of  a  Hermite- 
Gaussian  function  with  a  nuclear  center  may  be  expressed: 

[f(n,a,A;£)  |  Z^]  =   (2TT/a)  Z^  R(n,  a ,  A  -  £)  (A-40) 

Using  the  relationship  expressed  in  Equation  A-23  we  may  also  express 
the  electrostatic  interaction  integral  of  the  product  of  two  Hermite- 
Gaussian  functions  with  a  nuclear  center: 
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[f(n,a,A;_r)  f  (n  '  ,b,B;  r_)  |Z  J  =  — il 1        z     g  (K,n  ,n' ,a,b) 

'-  (a  .  b)   i  ^ 


A    A       A 


f(K,ab/(a+b),A-B^;£)  R(n+n '-K,a+b,Q-C_)  (A-41) 

where  ^  is  the  vector  defined  in  Equation  A-8. 

A-^.  Electron  Repulsion  Integrals 

The  electrostatic  interaction  integral  of  two  Hermi te-Gaussian  func- 
tions with  each  other  may  be  expressed: 

3„  r.3.,  1 


[f(n,a,A;r_)  |  f  (n",c,£;£' )  ]  =  /d^r  /d\ 


r-r ' 


f(n,a,A;r_)  f  (n",c,C_;£' )  {k-hl) 

We  may  again  use  the  relationships  derived  by  Boys  (1950)  for  spherically 
symmetric  Gaussians: 

[f(0,a,A;£)  I  f(0,c,£;£')]  =  A  F^ (  acT^/[a+c])  (A-43) 

where  X  is  now  the  new  constant 

\   =   2tt^^^  /  {ac  (a  +  c)'^^}  (A-A^t) 

and  the  vector  J_   is  defined: 

T  =  A  -  C_  (A-45) 

The  methods  used  in  Section  A-3  for  evaluating  the  nuclear  interaction 
integrals  may  be  extended  trivially  to  yield: 
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[f(n,  a,  A;  r)  |  f(n",  c.  £;  r ' )  ]  =  ("0    {2Tr^^V[ac  (a+c)"^]} 

R(n+n",  ac/(a+c),  A  -  £)  (A-46) 

Similarly  Equation  A-23  may  be  used  to  evaluate  the  electrostatic  inter- 
action integral  of  two  Hermi te-Gaussian  functions  with  a  single  Hermite- 
Gaussian  function: 


[f(n,  a,  A;  r_)f(n',  b,  B_;  _r)  |  f(n",  c,  £;  r_')] 

=    (-1)""  {2Tr5^^/[(a+b)c  (a+b+c)^^^]}  I   g(K,  n,  n',  a,  b) 


f(K,  ab/(a+b),  A-B_;  O)    R(n+n'+n"-K,  (a+b)c  ,  Q-£)    (A-47) 

a+b+c 

This  procedure  could  be  extended  straightforwardly  to  the  case  of  two 
Hermi te-Gauss ian  functions  interacting  with  the  product  of  two  other 
Hermite-Gauss ian  functions.   This  has  not  been  done  since  no  integrals 
of  this  form  are  necessary  for  the  current  work. 


APPENDIX  B 
MULT  I  POLE  EXPANSIONS  OF  COULOMB  INTEGRALS 


B-1.  Multipole  Expansions  of  Electrostatic  Interactions 

Consider  the  sum  of  electrostatic  interaction  terms  of  a  charge  den- 
sity P,(£.)  centered  about  the  origin  and  a  charge  density  P,(r_)  centered 
about  a  lattice  site  R^  of  a  two-d imensional 1 y  periodic  Bravais  lattice: 

U  =  I   /d^r,  Jd^r^  p,(l,)  02(12)       1  (B-1) 

where  the  coordinates  r_.  and  r_     are  illustrated  in  Figure  B-1.   These 
charge  densities  are  assumed  to  be  localized  about  their  respective  cen- 
ters, such  that  p  ( r_)  goes  to  zero  faster  than  any  polynomial  of  |r_|  in 
the  asymptotic  limit  as  |r_|  tends  toward  infinity.   This  condition  may 
be  more  restrictive  than  necessary,  but  it  is  easily  satisfied  by  the 
Gaussian  functions  used  in  this  work  as  well  as  by  exponential  functions. 
If  this  expansion  is  summed  in  order  of  increasing  |r|,  the  convergence 
of  the  Coulomb  terms  is  usually  slow  so  that  many  terms  must  be  included 
in  a  truncated  expansion  so  that  the  sum  may  be  computed  within  a  spec- 
ified precision. 

Piela  and  Delhalle  (1978)  have  found  that  the  use  of  a  multipole 
expansion  as  an  approximation  for  all  terms  having  [R_|  greater  than  some 
radius  R'  for  one-d imens ional 1 y  periodic  systems  such  as  polymers  is  an 
effective  method  for  speeding  the  computational  procedure.   The  multi- 
pole  expansion  considered  in  this  work  uses  the  bipolar  expansion  of  the 
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Figure  B-1.  Coordinate  system  for  the  multipole  expansion. 


inverse  of  the  interelectron  ic  separation  \r_.    "  £o  -  R_I   ,  as  presented 
by  Steinborn  and  Ruedenberg  (1973): 


— /   —     L,tn  K.   ,Tn 

(B-2) 


The  functions  ^(O  and  2 /> (_r )  are  the  modified  regular  and  irregular 
solid  spherical  harmonics  respectively  (all  notation  and  phase  conven- 
tions are  consistent  with  those  used  in  Steinborn  and  Ruedenberg,  1973) 

i^(r)    =    r^  y|^(9.<},)    [(2£+1)    (£-m)  !    (£+m)  !    /   k-n]~^^^  (B-3) 

and 

I'liL)   =   r"^^""^    Yj(e.<^)    [iiTT    {l-m)\    {t+m)\    /    {ll+\)]^^^  (B-M 


Equation  B-2  will  be  valid  for  values  of  r_,  and  _r„  for  which  one  of  the 
following  two  conditions  is  satisfied: 

')   111  I  <  IL2  ■*"  il  3"^  IL2'  '^    l-l'  °'' 

2)   IL2I  <  111  -  il  3nd  I^J  <  |R|. 
The  multipole  expansion  technique  approximates  the  electrostatic  inter- 
action with  the  truncated  expression: 

/d^r,  /d^r.  Pi  (ill)  ^2(1-2^   =   Z'""^  l   l^iR)      E        E 

IL,  -  L2  "  -'    '-"°    '^^      -^"^   m=max[-£,M+£-L] 

(-D^+W--^  M,(£,m)  M2(L-£,M-m)  (B-5) 


where  the  expressions  M. (£,m)  and  M_(L-£,M-m)  refer  to  the  multipole 
moments  of  p,  (r_)  and  p.,  (£_)  : 

M.(£,m)  =  /d^  i^(0  p.(r)  (B-6) 
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Stolarczyk  and  Piela   (1979)  have  demonstrated  that  the  truncated  expan- 
sion in  Equation  B-5  will  be  invariant  to  translations  and  rotations  of 
the  coordinate  system  as  long  as  |r|  remains  constant.   Since  all  omitted 
terms  in  Equation  B-5  for  the  multipole  expansion  will  be  of  the  order  of 
|R|    max    in  the  asymptotic  limit  as  |R_|  tends  toward  infinity,  the 
convergence  of  a  lattice  sum  using  the  truncated  expansion  implies  that 
the  exact  expansion  will  also  be  convergent. 

Thus  we  will  choose  an  approximation,  to  any  desired  precision,  of 
the  exact  lattice  sum  U  in  Equation  B-1.  This  approximation  v-yill  be  of 
the  form: 

u'  =    S   jd\.    jd\.    Pl^ll^  P2^il2^ 
|r1<R'  ^  |_r,  -  12  "  il 

+     E    y!^^""     E  if(R)  I     E  (-l)^""'^"-^  M,(£,m)M_(L-£,M-m) 
|R_|>R'   L=0   M   '-         1=0  m  '      ^ 

(B-7) 

A 

where  R'  is  chosen  according  to  some  as  yet  unspecified  criterion.   This 
leads  to  a  computational  savings  since  this  sum  may  be  evaluated  as  two 
terms:  One  term  is  the  sum  of  multipole  expansions  terms  for  all  lattice 
sites  exce pt  the  origin  (R_=  0_): 

U^  ,,  =  E   E{   E  i^(R)}  E   E  (-l)'-'^''"-^  M,(£,m)M„(L-£,M-m)  (B-8) 
™'^    L     M   |R|>0  '--      I     m  '      2 

vyhere  the  indices  range  over  the  same  intervals  as  those  given  in  Equation 
B-5.   The  value  of  R'  is  chosen  by  summing  in  order  of  increasing  |R_]  the 
difference  between  the  exact  electrostatic  integral  contained  in  the 
left-hand  side  of  Equation  B-5  and  the  truncated  multipole  expansion  In 
the  right-hand  side  of  Equation  B-5  until  the  difference  is  less  than 
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some  specified  tolerance.   This  sum,  including  the  exact  electrostatic 
integral  for  the  zero  lattice  vector,  yields  the  expression: 


'diff 

1-1   -2 


r.  -  r„  I        IR.|>0 


{  jd\.    /d^r.  Pl^ill)  P2(J:2^   -   E   Z  i';'(R) 

|L]  -  L2  ■  i|     ^  '^' 
E  E   (-l)'-"*'^'"^  M,  (£,m)  M-(L-£,M-m)}         (B-9) 


Summing  the  expressions  in  Equations  B-8  and  B-9  results  in  the  expression 

given  in  Equation  B-7.   The  evaluation  of  U   ,   wi I  1  be  considered  shortly. 

mu  1 1 

B-2.  Evaluation  of  the  Modified  Irregular 
Solid  Spherical  Harmonics 

The  lattice  vectors  R_will  have  a  zero  component  in  the  z  direction 
since  the  z  direction  is  chosen  as  the  normal  to  the  plane  of  periodicity. 

''171 

Thus  the  modified  irregular  solid  spherical  harmonics  2 « ( R_)  may  be  ex- 
pressed : 

i"l{R)   =  (-1)"  U-m)\    P^(0)  r'^^-"')  e'""*  (B-10) 

where  if  defines  the  angle  between  R^  and  the  x  axis.   The  Legendre  poly- 
nomial P/7(0)  will  be  nonzero  only  for  Z-m   even,  such  that  t-m  =   2k.    where 
k    is  an  integer.   if  this  condition  is  satisfied,  the  Legendre  polynomial 
wi  1 1  equal : 

P£(0)  =  (-1)^^  {l+m)l    /   [2^  k\    (fe+m)!}  (B-ll) 
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Since  every  Bravais  lattice  will  possess  inversion  symmetry,  we  may  con- 
sider tlie  sum  of  modified  irregular  solid  spherical  harmonics  for  both 
R  and  -R^.   Thus  we  may  utilize  the  relationship 

to  yield  the  expression 

i"l{R)   +  Z^-R)  =  (-;)^  (^^^)'  ^^-11     e""*  (B-13) 

^~  ^     ~  2^"'  fe!  (fe+m)!  R^  ' 

for  even  values  of  t   and  m  such  that  2k   =  £-m.   if  either  £  or  m  is  odd 

then  this  sum  instead  will  equal  zero.   In  the  present  work,  L  has  the 

^  max 

value  of  six  and  thus  fifteen  real  numbers  must  be  evaluated  for  the 
values  of  the  modified  irregular  solid  spherical  harmonics.   These  fif- 

teen  real  numbers  correspond  to  the  real  values  of  2_,  Zi  »  ^""^  ^c,'    ^""^ 

"2      ^2      "k      "2      "U  ^6 

to  the  real  and  imaginary  parts  of  Z- ,  Z , ,  Z, ,  Z^,  Z,,  and  Z,.      No  value 

is  computed  for  Z^  since  no  multipole  expansions  will  contain  a  nonzero 

monopole-monopole  interaction.   No  additional  values  are  computed  for 

Ip      since  these  may  be  generated  from  the  relationship 

Z;^'"(R)  =  (-1)''[Z'J(R)]"  (B-14) 

The  value  of  the  sum  of  ^p^^    ^°''  ^^^  1^1^^  '^  computed 

by  summing  all  lattice  site  contributions  for  a  finite  number  of  lattice 
sites  up  to  a  maximum  radius  R".   The  sum  of  lattice  site  terms  for  |R_| 
greater  than  R"  is  approximated  by  the  integral: 
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,  ,^   Z/(R)  ~-   ^'^    /:..  r  dr  /2-  i;(r)  d*  (B-15) 

=  (-l)'^2ir(£!/fe!)^5  „  /  {2^nU-^){R"^~^)} 

mo 

where  fe=£-m  and  Q    is  the  area  of  the  two-dimensional  unit  cell  of  the 
two-dimensional  Bravais  lattice. 


B-3-  Evaluation  of  Multipole  Moments 

The  evaluation  of  the  multipole  moments  of  the  nuclear  charge  is 
uncomplicated.   A  nuclear  center  having  a  total  charge  Z.  located  at  a 
position  vector  A^will  have  the  multipole  moment  M(£,m): 

M(£,m)  =  Z.  P^(A)  (B-16) 

The  multipole  moments  of  a  charge  distribution  represented  by  a  Hermite- 
Gaussian  function  f(n,  a,  A^;  r)    may  be  evaluated  using  the  relationship 
for  modified  regular  solid  spherical  harmonics  (Steinborn  and  Reudenberg, 
1973): 

t 


Setting  r_.    =  £_  "  A^  and  £_„  =  A^,  we  see  that 

jd\    f(n,a,A;r)  y^^(r)  =    ^^   ^      ^      ^£-r  (^)  1^^'    ^i  (j:-A)exp  [-a  |  r^-A  |  ^1 

—  t'   m' 

(B-18) 

Since  the  Gaussian  function  is  spherically  symmetric  about  A^,  the  inte- 
grals on  the  right-hand  side  of  Equation  B-18  will  be  nonzero  only  for 
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the  term  where  £'  is  equal  to  zero.   Using  the  definition  of  the  modified 
regular  solid  spherical  harmonics  given  in  Equation  B-3,  we  may  evaluate 
the  integral  for  the  term  where  Z'    is  equal  to  zero: 

jd\   f(n,a,A;r_)  /^^(O  =  9^  V^(A)  U/a]^^^  (B-19) 

If  the  charge  density  p  ( r_)  has  the  form  of  the  product  of  two  Hermite- 
Gaussian  functions: 

p(r_)  =  f(n,  a,  A;  r_)  f{n',  b,  B_;  r)  (B-20) 

we  may  use  the  relationship  expressed  in  Appendix  A  in  Equation  A-23  in 
conjunction  with  Equation  B-19  to  yield  the  multipole  moments  of  charge 
densities  which  have  the  form  described  in  Equation  B-20: 

M(£,m)  =  E   (K,n,n',a,b)  f  (K,ab/[a+b]  ,A-B_;0_)  [Tr/(a+b)]^' 
K 

^p'"^^^(Q)  (B-21) 

where  Q  is  as  defined  in  Equation  A-8. 

The  functions  '  o  ( A_)  may  be  evaluated  by  considering  the  following 

m  '  £ 

expansion  of  the  spherical  harmonics  Y„(r)  multiplied  by  a  factor  of  r  , 

where  m    is  greater  than  or  equal  to  zero: 

I  ymr.    _    {-]f/   r(2£+l)    (£-m)'    1/2      'uH    n      2.m/2  d^"""      ,2    ,.£ 
^  2^  II  kv        U+m)  !  d^^  '" 

(B-22) 

where  r,   =   cosO.   We  may  expand  a  factor  of  the  above  equation  in  the 
fol lowing  form: 

r  (1-?  )    e   ^  =   X   [,J  1   y  X  (B-23) 

fe=0 
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We  may  also  expand  the  factor: 

d?  O^lj^l-m  ^    [l-lj-m) ! 

Thus  by  combining  the  expressions  in  Equations  B-22  through  %-lk   and 
using  the  definition  of  the  modified  regular  solid  spherical  harmonics 
given  in  Equation  B-3,  we  find  the  following  expansion: 

^      2^£!(£+m)l   j     k       ^       i      {l-li-m)\ 

r  -'  X    y  z   -^  (B-25) 

We  may  also  expand  r  -^  into  products  of  powers  of  the  Cartesian  coor- 
dinates X,  y,  and  z: 

r^^'  =  {  \        [i]     f  1  x2^  y2P-2^  z2^2p  (^.26) 

Thus  we  have  an  expansion  of  '»(£_)  in  products  of  powers  of  the  Cartesian 
coordinates: 

^(L)  =  t'^^ ,    ^  ^  ^  M^)  (f)  g  [p  i' (-!)'■ 

^     2^£!  U+m)  !   j  p  q  fe      ^   '^  ^ 

{2t-2i)\  ^m-h+lq   ^2p+k-2q   _l-lp-m  ^^_^^^ 

(£-2y-m) ! 

We  may  trivially  evaluate  the  quantity: 


101 


^A  ^(A)  =  /-;)'     z  z  z  zi^(-i)^'  Q  6)  0  0 

-       ^  2^11  U+m)  I      j     p     q     k  t^       J        P       Ci 

(2£-2/) !  {m-k+2q) !  {2p+k-2q) I  (£-2p-m) ! 

U-2j-m)\      (m-fe+2Q-ni)!      (2p+fe-2q-n2) !      (£-2p-m-n3)! 

^m-fe+2q-ni      y2p+fe-2ci-n2     ^£-2p-m-n3  ^g_2g) 

where  the  summation  indices  have  the  following  ranges: 

0  <  2j   <  I  -  m  (B-29) 

£  -  m  -  n3  -  2p  -  nj  +  n2  -  m                            (B-30) 

2p  -  n2  +  m  ^  2q  <  ni  -  m  (B-3I) 

2q  +  m  -  ni  ^  fe  <  2q  +  n2  -  2p                           (B-32) 

This  set  of  equations  completes  the  evaluation  of  all  necessary  quanti- 
ties for  the  computation  of  the  multipole  moments  of  Hermi te-Gauss Ian 
functions  and  of  nuclear  centers. 


APPENDIX  C 
RECIPROCAL  LATTICE  EXPANSIONS 


C-1.  General  Properties 

Many  of  the  quantities  evaluated  in  this  work  may  be  expressed  as 
an  expansion  in  the  lattice  vectors  Regenerated  from  the  primitive  lat- 
tice vectors  R^,  and  R^  defined  in  Section  2-1: 

f(£)  =  S  exp[il^  -R]  U(£+  R)  (C-1) 

R       " 

where  _r  i s  a  vector  in  a  three-dimensional  Cartesian  coordinate  system. 
Let  us  define  £  then  by  the  Cartesian  coordinates  r_=    (x,  y,  2)  and 
furthermore  define  s^  to  be  a  vector  in  a  two-dimensional  Cartesian  co- 
ordinate system  such  that  s_  =  (x,  y) .   Since  the  Bloch  function  f ( 0 
has  the  property: 

f(£+  R)  =  exp[-ik//-R]  f(_r)  (C-2) 

we  may  express  f (0  as  the  product  of  a  phase  factor  and  a  periodic 
function  f  (s_,  z)  such  that: 

f(0  =  exp[-ik//-s]  f  (s_,  z)  (C-3) 

where  f  (s_,  z)  is  periodic  only  in  the  coordinate  s_: 

^'p^l^  R.  z)  =  fp(l,  z)  (C-i*) 
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We  may  expand  f  U,  z)  as  an  expansion  in  the  reciprocal  lattice  vectors 
K  described  in  Chapter  2  as: 


f  (s_,  z)  =  E  exp[iK-s_]  f(K,  z)  (C-5) 

P         K 


where  f(j<,  z)  may  be  expressed  as  an  integral  over  the  entire  plane 
defined  by  s^: 

f(K,  z)  =  fi"^/d^s  exp[-i(K-k^^)-£]  [){r)  (C-6) 

Thus  we  may  expand  our  original  function  f (£)  in  a  reciprocal  lattice 
expansion  of  the  form: 


f(£)  =  Z  f(K,  z)  exp[i(K-k   ).s_]  (C-7) 

K  ^^ 


As  an  example  of  this  approach  let  us  consider  the  reciprocal  lat- 
tice expansion  of  a  Bloch  function  constructed  from  the  Hermi te-Gauss ian 
functions  f(n,  a,  A;  r)    defined  in  Appendix  A: 

*(D  ii//)  =  ^  exp[ij<^^-R^]  f(n,  a,  A;  £  +  R_)  (C-8) 

R 

The  expansion  coefficient  f  (j< ,  z)  will  then  equal: 

f(K,  z)  =  f^'Vd^s  exp(-i(K-l<^^)-s_]  f(n,  a,  A;  r)  (C-9) 


^  3^  exp[-i(K-l<   ).A]  exp[-a(z-A  )  ] 


Tr_  j-Cni+n,)/,,  _.    \^\(u      ,    Nn,  no/Z 


HpgCa^  ^(z-A^)]  exp[-a(z-A^)^]  exp  [- i (K-k^^) -A] 
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Thus  we   see   that 


J./        I       ^        ■'^      .-(ni+no)      n3/2  r    1/2,      .    v,  r      /      ,    n?i 

(filr;    k.,)    =-75-1        ^      ^'    a    ^        H      [a        (z-A   )]    exp[-a(z-A   )^] 
—     — /  /ail  n  3  z  z 


K 


Overlap  matrix  elements  may  be  calculated  using  Equation  C-10  if  we  use 
the  identity  (where  the  integration  over  £  is  over  the  entire  three- 
dimensional  space): 

Z  exp[ik,,-R^]  /d^r  f  (n,  a,  A;  r)    f{n;  b,  B^;  £  +  R^) 
R       ^^ 


=  [^]^^^-l)"  Z  exp[ik^/-R]  f(n  +  n';  ^,  B-A;  R)       (C-ll) 
K 


since  this  expression  is  equivalent  to  a  Bloch  function  i)  ( r_;  k^//)  as 
defined  in  Equation  C-8  evaluated  at  r  =  0,  multiplied  by  a  constant. 


C-2.  Coulomb  Integrals 

Much  of  the  work  using  reciprocal  lattice  expansions  to  expand 
Coulomb  integrals  relies  upon  the  expansion  of  the  Coulomb  potential  in 
a  reciprocal  lattice  expansion.   One  difficulty  with  this  approach  is 
that  the  long  range  behavior  of  the  Coulomb  potential  leads  to  many 
problems  with  the  K_  =  £  term  in  the  reciprocal  lattice  expansion,  although 
Schwalm  and  Monkhorst  (I98O)  have  pointed  out  how  one  may  treat  the 
difficulties  in  this  approach.   A  more  useful  approach  for  the  purposes 
of  this  work  is  first  to  create  a  direct  lattice  expansion  of  the 
desired  integral  and  integrate  each  term  analytically.   We  now  proceed 
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to  construct  such  an  expansion  and  derive  an  equivalent  expansion  in 
reciprocal  lattice  vectors. 

First  let  us  consider  a  Coulomb  integral  of  the  form: 

V(s_)  =  E  {/d^r  Jd^r'  exp  [-a  |_r-£-A-R  |  ]  exp[-b  |_r ' -B^|  ]  ,  |^,  i 
R  —  — 

N  1 

-[^]^^^  ^     E  Z^  /d^r  exp[-a|r_-s_-A-R|^]  i  _^  i  }    (C-12) 
tot  m=  1  ' m' 


where 


N 

tot      ,   m 
m=l 


It  has  been  demonstrated  in  Appendix  B  that  this  lattice  sum  will  be 

convergent  if  we  sum  in  order  of  increasing  |  R_|  .  We  may  reduce  Equation 
C-12  to  the  form: 

V(s)  =   T.    {      I      a      j   exp[-3  I  s+r  +R|^u  ]du}  (C-l4) 

where  terms  for  which  the  index  m  is  zero  correspond  to  electronic  terms 
such  that: 

ao  =  2u^^V{ab[a+b]^'^%  (C-15) 

3o  =  ab/[a+b]  (C-16) 


lo  =  A  -  B  (C-17) 
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Terms  for  which  the  index  m  is  greater  than  zero  correspond  to  nuclear 
terms  such  that: 

a     =   -  277^/^  Z  /{Z^  ^  a  b^''^}  (C-l8) 

m  m   tot 


^m  =  ^  (C-19) 


r  =  A  -  C  (C-20) 

— m   —   — m 

We  may  expand  VCs)  in  a  reciprocal  lattice  expansion  of  the  form: 


V(s_)  =  Z  exp[iK-s_]  V(K)  (C-21) 

K 


where  each  coefficient  V  (j<)  is  given  by: 

V(K)    =   Q~      /d^s   exp[-iK-£]    V U )  (C-22) 


N 

=  Q        /d   s   e    '-'—   Ldu      E      ex.  exp[-6    Is+r    I    u   ] 

•'  ■'0      .   I       m' m' 

m=0 


For  lattice  vectors  such  that  \k\    is  greater  than  zero,  we  may  inter- 
change the  order  of  integration  in  Equation  C-22  without  difficulty. 
Bearing  in  mind  that  restriction,  we  interchange  the  order  of  integration 
to  yield: 

V(K)  =  TTfi"^   E  -^exp[iK-r  ]/^u"^exp[-6  z^u^-kV'^S  u^]du    (C-23) 

—  „  6    ^ m  •'o      "^   m  m        m 

m=0  m 

A  proper  treatment  of  V(j<)  requires  a  more  complicated  analysis  for  the 
case  of  ]<  =  £.   The  reason  that  the  order  of  integration  may  not  immedi- 
ately be  interchanged  is  the  singularity  of  the  integral  over  s_  for  values 
of  the  integrand  near  u  =  0.   By  use  of  the  expression: 
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/     exp[-6    |s+r    pu^]    du 
0  m m 


zF^—J      ""z^"^     ^iirrr    s^ 

m  m  ' m ' 


-/: 


exp[-6    |s+r    I      u    ]    du 


{C-2h) 


let  us  examine  the  integral 


-1  r.2  ,,y       1         "m        1/2, 


m=0  R^'^  m=0  '^   ^m     '-  -i-ml 


2   2, 


-  (  Z  a   /  du  exp[-B  |s+r  |   u  ])} 


(C-25) 


m=(. 


-1/2 
term  by  term.   We  see  that  the  values  of  the  quantities  a  6      for  the 

m  m 

case  of  m  =  0  yield 


„  -1/2   .  5/2   ,n-3/2 
ao6o     -   ^-rr    LabJ 


(C-26) 


and  for  the  case  of  m  greater  than  zero 


•1/2 


Zir   i. 


m 


mm 


^otf^^^ 


3/2 


(C-27) 


lead  to  the  result: 


'   «mC'^'  =  ° 
mm 

m=0 


(C-28) 


in  order  to  evaluate  the  contribution  of  the  second  term  of  Equation  C-25 
let  us  consider  the  integral: 


/ 


^2    r  1 

' m ' 


1  1   r^2  ,  1 


1 


•} 


s+z  k 
'—  m  ' 


(C-29) 


=  U  fC 


0  [s^z  ']'^' 
m 


1)  d! 
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We  may  transform  this  integral  to  the  form: 


^-r(, .  Li/2-  ^)  '^- 


0        [s^+2     ^] 

m 


-     2TfZ        / 

m   •'o 


Tr/2  de 


1    +   sin( 


(C-30) 


By  use  of  the  indefinite  integral 


r   de 

•'  1  +  sine 


tan  (-jj-  -  j) 


(C-31) 


Equation  C-29  reduces  to 


/«'=  (Til^  -  4r' 


27rz 


m 


I ml 


(C-32) 


Interchanging  the  order  of  integration  of  the  last  term  of  Equation  C-25 
leads  to 


.2   f" 


N 


2   2, 


fd  s  /  du   E  a  exp[-B  js+r  ]   u  ] 

•'     ■'i       „   m       m m 

^     m=0 


m  (•<»  r  „   2  2,  du 

=  TT  Z  ^-  /  exp[-6  z  u  J  — 

„  D   1  mm     2 

m=0  m  u 


(C-33) 


Integration  by  parts  yields 


f°°  r    n      2    2t    du  r    o      2-1 

J      exp[-6   z   u    J   -—  =   expl-e   z    J 
•'i                    m   m             2  mm 

^  u 


2    r°° 


m  m 


2B  z      J     exp[-6  z  u-"]    du      (C-3A) 
m   m   -^l  mm 


By   use  of    the   relationship 


/;  exp[-B^z^u^]du   =  ^i-^)'''    -    (   exp[-B^z^u2]du 


(C-35) 


i    Z 

m  m 


-    W      ^      ^^/-    -    F    (ft   z^) 

6  z 
m  m 
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where  F^  Fs  the  incomplete  gamma  function  defined  in  Appendix  A,  we  may 
reduce  Equation  C-33  to  the  form 

/d  s  /;du  J  a^   exp[-6js.r^|2  u']    =  .   z     a  {-L  exp[-B  z'] 


m=0  '"       I"  --m-  m^B 


m=0 


m 


m  m 


/  7T  »  1/2  o  , 

-  {-^J  z   +  2z   F  fB  z  )}  /r  ,^n 

Pn,      m     m  0^  mm"^-^  (C-36) 

Notice  that  while  the  integral  in  Equation  C-35  will  be  undefined  for 
z^  =  0,  an  analysis  of  Equation  C-33  indicates  that  Equation  C-36  will 
remain  valid.   Inserting  the  results  of  Equations  C-29,  C-32,  and  C-36 
into  Equation  C-25  yields 


'  -^  ^  "  '""'   '   "m^-^  f  ^^P^-^yj    ^  2zV„(B  z')}  (c-37) 

m=n    '^    p_        mm       m  0  ^  m  m' -^  VL  i/) 


Let  us  consider  the  evaluation  of  the  followin 


g  expression: 


^A  i  V(l)  =  I   e'i<-l  9^  9^'  v(K) 


(C-38) 


The  partial  derivatives  with  respect  to  the  vector  components  B.  do  not 
apply  by  construction  to  nuclear  terms.  Since  V ( 0 )  depends  only'upon 
'^x'  \'  ^x'  °''  \  ^^  ^^^  that 


nD  ^n ' 


^  i    -      "  nj.O  no  ""n-.O  ^n'  0 


[(-,)"3  ,^(__,(,^,^,^  3^^  ^^.  ^)  ^  2^2  ^(^^^..^  ^^^  ^^^^ 


N 

"3'0  i,  "-^t  '^"z'  ^'  ^n.=  i)  -   K   ^(^'  B^'  O)]  (C-39) 


no 


where  the  notation  n  denotes  the  set  (O,  0,  n  ).  For  vectors  J<  not  equal 


to  the  zero  vector: 


6"  5"'  \/(K)  =  ^  {il  i^l+na-nl-n^  ^  n^-n^  ^   n2-n'  ^  ng+n^ 
^  i    -    2  Bg  X      y       Zq 


{/     exp[-6,z2u2-K2AB„u2]   ^}  +  5    ,    .    5    ,    ^   6    , 

•^0  0    0  0  ^2  "1.0      n2»°      "3' 

N      a  1 

m  .ni+n9    l'   "^1    1/   "2    '^^3    f  r    o      2    2    i/2//.r>      2i    du   , 

E     -—  I     ^      ^    K      '    K     ^8^1      expl-B   z^u'^-K'^/hB   u^J   —  } 
,B  xyz-'omm  m  o' 

m=l      m  m    "  u'^ 


(C-40) 

We  may  reduce  the  integrals  in  Equation  C-'40  to  known  form  if  we  consider 
the  following  construction. 

First  let  us  define  the  quantity  fo(a,b): 

fo(a,b)  =  /^  exp[-a2u2  +  b2/u2]  ^  (C-4l) 

0  ^2 

=  -^TT-  (e     /  (a  +  — )  exp[-(au-b/u)2]  du 
2b       '0  ^2 

-  e    /  (a  -  — )  exp[- (au+b/u)  •  ]  du  } 
0      ^2 

1  r  ~2ab  (■■»    -V-'  ,   ,   2ab  /•«>    -v2 

=  2b  ^^     Jb-a  ^   ^^  ""  ^    ■'b+a  ^   ^^  ^ 


Since 


r  e'^'dv  =  ^'^'-  f  e'^'dv  (C-'42) 


TT 


1/2 


=  T   "  y'^o(y^) 


2 


11 1 


we   see    that 


fn(a,b)    =  j^  {TTl/2cosh   2ab   -   e^^^(b+a)    FQ([b+a]2) 


e"^^^b-a)    Fo([b-a]2)} 


(C-i^B) 


Consider  the  relationship: 


J-   (e^^^^bia)  /^  exp[-(b±a)2u2]du}  =  i   {2b  e^^^^bia)  Fq  ( [b*a]  2) 


3a 


+  e*2^^   Fo([b±a]2)    -    2e*2"^b±a)2   /'u2e- ^^^^^ '^'du}  (C-H) 


Utilization   of    the    identity 


A    2    "Y^u2               1     rA      "Y^u^  -y    , 

J    u'^e  du  =  Lj      e  du   -   e        J 


2y 


2      0 


(C-^5) 


results    in 


^  {e*2^^b±a)    /     exp[-(b±a)2u2]du} 
da  0 


=  ±  2b{e-2^'^(bia)    Fo([b±a]2)}  *  e"^^'^^'^ 


{C-k6) 


Thus  we  see  that 


fl(a,b)    =  -^  fo(a,b)    =   7rV2sjnh   2ab    -   e^^^b+a)    Fo([b+a]2) 


+  e"^^^b-a)    Fo([b-a]2) 


ic-kj) 


and 


f2(a.b)    =  4b2    fo(a,b)    -   2e"^^'"'^'^ 


(C-i48) 
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The  above  results  then  lead  to  the  recursion  relations 


f  (a,b)  =  a"  fo(a,b)  ic-hS) 

n         3 


n-2 

^a 


n  -c  ,.  /   ,  \ 
=  9,    ^2^^'^^ 


=  ^b2  3-2  f  (a,b)  -  29"-2  e-^^'-^^') 

3      U  3 


=  I.b2  f   (a.b)  -  2(-l)"  H  ,(a)  e"^^'"''^'^ 
n-2  n-2 

where  H  (x)  refers  to  the  Hermite  polynomials.   This  reduces  the  factors 
n 

in  Equation  C-'tO  to  the  form: 


9"   /   exp[-6  z2u2  -k2As  u2]  ^ 
z^  0        mm         "1    u2 


m      n   m      m    '  — '     m 


m    n  m    m 


APPENDIX  D 
NUMERICAL  INTEGRATION  GRID 


We  wish  to  describe  briefly  the  scheme  for  selecting  the  numerical 
integration  grid  points  and  weights  used  in  the  exchange  fitting  proce- 
dures discussed  in  Chapter  III.   A  general  review  of  numerical  integration 
techniques  is  contained  in  Davis  and  Rabinowitz  (1975)-   Several  dif- 
ferent calculations  have  been  performed  testing  various  schemes  for  numer- 
ical integration  over  the  three-dimensional  unit  cell  defined  in  Chapter 
II,  with  results  indicating  that  the  details  of  the  selection  scheme  are 
not  as  important  as  the  total  number  of  points  chosen.   Two  basic  guide- 
lines have  been  followed  in  all  schemes  tested: 

1)  The  arrangement  of  grid  points  near  each  nuclear  center  should 
remain  the  same  with  respect  to  changes  in  the  lattice  constant. 
This  condition  minimizes  the  potential  problems  arising  from 
changes  in  the  grid  in  regions  of  high  density  near  the  nuclei, 
which  would  tend  to  mask  the  effect  of  actual  changes  of  the 
density  on  the  desired  integrals. 

2)  The  coordinates  for  the  three-dimensional  grid  should  be  chosen 
so  as  to  take  advantage  of  the  property  that  Hermi te-Gauss ian 
functions  may  be  factored  into  the  product  of  independent  func- 
tions of  each  of  the  Cartesian  coordinates.   This  property 
results  in  being  able  to  factor  the  Bloch  functions  which  we 
must  calculate  at  each  grid  point  into  the  product  of  a  function 
of  the  coordinates  parallel  to  the  surface  (x  and  y)  and  a 
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function  of  the  coordinate  normal  to  the  surface  (z) .   Choosing 
a  grid  where  several  points  share  the  same  value  of  the  coor- 
dinate normal  to  the  surface  for  each  choice  of  the  normal  co- 
ordinate can  lead  to  a  substantial  savings  in  computation  time 
in  evaluating  the  Bloch  functions. 
The  scheme  used  for  the  calculations  presented  in  this  dissertation 
begins  by  constructing  a  grid  using  the  cylindrical  coordinates  (p,4),z) 
about  each  nucleus  located  in  the  primitive  unit  cell  under  consideration. 
Each  grid  point  lies  at  the  center  of  a  representative  volume,  with  the 
weight  assigned  to  that  sampling  point  equal  to  approximately  that  part 
of  the  representative  volume  contained  in  the  region  of  space  nearer  to 
the  central  nucleus  than  to  any  other  nucleus  (including  nuclei  in  ad- 
jacent unit  cells).   The  values  of  the  normal  coordinate  (z)  were  chosen 
using  the  radial  distance  of  every  fifth  point  of  the  Herman-Ski  1 Iman 
mesh  (Herman  and  Ski  1 Iman,  1963).   This  Herman-Ski  1 Iman  mesh  contains 
k^]    radial  points  where  the  interval  between  adjacent  points  doubles  at 
every  '♦O-th  point.   The  initial  interval  between  the  origin  (the  first 

point)  and  the  second  point  of  this  Herman-Ski  1 Iman  mesh  equals 

2      1/3 
0.0025(9iT  /128z)    ,  where  z  is  the  atomic  number  of  the  central  nucleus. 

The  values  of  the  radial  coordinate  P  were  also  chosen  using  every  fifth 
point  of  the  Herman-Ski  1 Iman  mesh  for  values  of  the  normal  coordinate  (z) 
sufficiently  small  in  magnitude.   At  progressively  larger  values  of  the 
magnitude  of  the  normal  coordinate  (z)  every  tenth,  fifteenth,  twentieth, 
etc.  point  of  the  Herman-Ski  1 Iman  mesh  was  used  for  the  value  of  the 
radial  coordinate  (p).   For  each  choice  of  the  normal  coordinate  and  the 
radial  coordinate,  four,  eight,  or  twelve  values  of  the  angular  coordi- 
nate ({i  were  used.   All  points  having  a  representative  volume  lying 
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completely  outside  the  region  nearer  the  central  nucleus  than  any  other 
nucleus  were  excluded  from  the  integration  grid  since  the  volume  assigned 
to  that  grid  point  either  lies  outside  the  primitive  unit  cell  integration 
region  or  will  be  treated  with  points  centered  about  another  nucleus  in 
the  primitive  unit  cell.   A  set  of  FORTRAN  subroutines  capable  of  calcu- 
lating the  sampling  points  and  weights  generated  by  the  above  mentioned 
procedure  is  given  in  Table  D-1. 
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Table  D-1.  IBM  System/370  FORTRAN  IV  code  for  generating  sampling  points 
and  weights  used  in  numerical  integration  procedures. 


SJTPQJT  INE     INTNUM(RhUCt      2CHG.      VtCT,      Nr  LC  ) 
IMPLICIT     PE^L* B( A-H, C-Z) 
LQ3IC^L     CJT,     TRUE.     ^ALSE 

c 

C  THIS     5U3PGUTI\E     CALCULATES     T  hE     MUMEPIC^L     GTI':     =20INTG 

C  ANO     WEIGHTS     USE-'J      IN     THE     EXCHASGc     FITTINj     pr-OCEDUnS. 

C  THE     SUIROUTINE     iNTGr^D     IS     REGUIPED     FOR     USE    OF     THIS 

C  SU-^ROUTI  ME. 

c 

C  IN'^UT     VARIARLES: 

C 

C       RNJC  (3t  NNUC  )  :     CARTESIAN  CCOROINATES  r:-  NUCLEI 

c 

C  ZCHG(NMUC):   ATC'IIC  NU^BSP  CF  NUCLEI 

C 

C  NM'IC:  NUV3ER  CF  NUCLEI   IN  PRI'^'ITIV'^  UNIT  CELL 

C 

C  yECT(2t2):        X     ANH    Y     CARTESIAN    COCRDINATES     OF     PP-r-iITIVE 

C  LATTICE     TRANSLATICNS 

c 

CI'IENSION     RNUC  (7  ,NNUC  )  t      VECTd.::),      ZChC(NNUC) 
OniNSION     INUC(IC).     ^(-^Z).     -^XC^Z),     WTZ(32) 
CATA     01/ 3.1  4159^  e535C575i5  +  00/ 

D^T\     T'-'UE  tFALSE/     .  T  i-U  E  .  t      .FALSe./ 
r>^TA      ^TOL/l.D-5/ 
CATA     NSK  IP/5/ 

■.yfRI  TE     (6,  ?0  JOO  ) 

IDT     =     0 

c 

C  Sn'^T     NUCLEI      IN     CROEP 

C 

no     20      1=1 f NNUC 
20  INUC ( I )      =1 

IF     (NNUC.EO.l)     GC     TJ     40 

DO    3  0     I =1 tNNUC 

OG     30     J=I .NNUC 

IF     {RrjJC(3.  INUC(  J)  )-PNUC(  ;.  INUC(  I  )  )  )      25.      EC.      20 
25  IN-     INUC ( I ) 

INUC (  I  )     =     INUC ( J  ) 

INUC(J)      -      IN 
30  CONTINUE 

C 
C  LCnP-CVfP     NUCLEAR     SITES      IN    i;nDE?J 

c 

40  on    960      IN=1,NNUC 

I  •'N     =     I  \UC(  IN) 

C  \LL      INTGKD(  /.CHG(  IZM  )  ,     NSKIP,     '<  .     RX.      '.>  TZ  .      NP\T.      ;C\Lr) 
'J^ll  T     =      :   .00      -      P(  2  ) 

c 

C  LOT-3     JVER     PAUI'VL     "CirJTS 


c 


on     9  50      l'^=l.NPNT 
OUT     =     TPU2 

")  E  L     =     P  {  I  P  )     -     R  (   I R  f  1   ) 
M3     =      4      V      I  O  INT  (ilEL/UNI  T  ) 
I^      (M)'.GT.12)      N-J     =      12 
IF     (^;".LT.4)     NP     =     4 
IF     (    r^  .(-0  .  1  )     NP      -     1 
\-7C     -    PI    /    .:flcat  (  rj=>  ) 
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Table    D-1 .     (Continued) 


C  LOO^     AROUND     C  I  HC  LV-FEFE  NCC     CP     Cir.'CLE 

C 

70  00     300      ic=l    ,NP 

ANGLE     =     ARC     *     0FLGAT{2*IP)      +■     PI      /     4  .  CO 

XO     =     =JX(IP)      *     CCCS(ANGLi) 

YO     =     RX(IS)      *      3S  IN(  A,NGL.:) 

XI     =      XO     +     PNUC  {  1  «  IZN  ) 

Yl      =     YO      -t-     RNUC  (J  .  IZN  ) 

WTO     =     1  . O  0 

c 

C  LCD=>    CVIZK     Z      (NCPVAL)      CI'^'CCTICN 

c 

NZ     =      -NPNT 

N2     =     2*MPNT     -      1 

00     3 CO     N=1,M2 

NZ     -     \Z     +     1 

ICHK     =     NSKIF     *      IflBS(NZ)      /     40      +     1 

ICHK      =      MINOflCHK ,     NPNT) 

IF     (MOD(IR-lt       ICH<).Nt.O)     Oi     TO     600 

ZO     =     OFLOAT(  I  S  I-.N(  1  ,NZ  )  )      *     t' x  (  I  ABS  (  N  2  )  f  1  ) 

Zl     =     ZO      +     RNUC  (  ~  .  IZM  ) 

CMX      =     0 .00 
IF     (    I^+  I^-MK.LT  .NFNT)      1^  M  X     =     ^  {  1 '^     *      ICHK) 

R>^X      =     -DLOG(      (R{IR)      +     RMX)/:-,DO)     /     SCAL'-I 

IF     (   IP-1  )      i:.'0«      1  jO,      140 
130  VGL=     PI      *     RMX**2 

GC     TO     150 
160  Rvi-g     =    -ULCG{      (RtlR)     +     R  (  I  P- I  CHK  )  ) /2  .  C  0  )  /  SC  AL  F. 

VOL     =     PI      *      (PMX      ♦     PMN)      *      ( RM X    -     RMN ) 
150  l»T     =     wTO     ^'     V.T7  {  I  AnS(  NZ  ) +1  ) 

VOL      =     VUL     /     DP LO AT (MP) 

V0L2     =     VCL 

RMX2     -     OSOPT(7  0**2      <-     PMX**2) 
C 

C  LOOP     OVER     AOJACZNT     UNIT     CULLS 

C 

DO     20D      IXX= 1 , 1 

n  T    200     I  YY=1  ,2 

XX     =     VtCT(  1  ,  1  )  *'-)hLQAT(  I  <X-2)      +     VEC  T  (  1  f  2  )  *^D  FL  C  AT  (  I  Y  Y-2  ) 

YY     =     VECT  (  2,  1  )  *DFLOAT(  IXX-2)      *■     VECT  (  2  t  2 )  ^OFLG  AT  (  I  Y  Y- 2  ) 
r 

C  LOOP     OVlP     NUCLEAR     C'^NT'£'<S      I\     UNIT     CELL 

C 

iOO     200     I  N2=l  ,NNUC 

IF     ({(rxX.EO.2)  .AND .(IYY.E0.2)  ).AND.   (IN2.EC.I2N>  ) 
X  GO     TO     200 

XI      =     "■(NUC(   1  .  I  i\2)      +■     XX     -     RNUCd.IzrJ) 

YI     =     PNJC(2,IN2)      +     YY     -     f^MUr(2,IZM) 

ZT      =     WNUC(:',IN2)      -     RNUC(3,rZK) 

^2     -     X  I  *  *  2     ■•-     V  !■>:  *  2 

R12     -     f:)SQRT('-i2     *-     ZI**2>     /     2.00 

IF     (  r  MX2.LT.R1  2)     GO     TC    .':  0  0 

IF     (  P2.  GT  .RTLL  )      CO     TO      1  r.  2 
C 
C  NUCLE-VR     CE^.TFr-     L  I '^  5      I  rj     GUPFACE     NCPVAL      CIPFCTICN 

c 

IF     ( 'JA3r;  (  ZO  )  .GT.CAGS  {  Z  I-ZO  )  )      GO     TO      ,-^0  0 
GO     TO     200 
165  FACTCn'     =      (1.00     -     Z  I  <- (  Z  I  -  2  .  OO  >>  ZO  ) /P  2  )      /     2.00 

v2     =     JGi.RT(R2) 
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Table   D-1 .    (Continued) 


C 

c 
c 


c 
c 
c 

17 


c 
c 

c 


175 

17c 

1  78 
200 

C 

c 

c 


c 

?00 

<;oo 
<;50 

9  6  0 


10  ) 
200 
200 


O    =     DABG(rACTGR     *     R2 ) 
AQUT     =      0.00 

IF  (D.GT.RMX)  GC  TO  200 
THZTA  -  OAPCCSO  /  ^-VX) 
THETP     -      O.OO 

IF     (D.LT.RN'N)      THETP     =CARCCS(D    /     R^N  ) 
CTSA     -     (XI*X0     +     ^fI<^YC)     /     (P2     *     WX(lr)) 
PS  I     =     0A9CCS  ( COS  a) 
IF     { FACTOP.LT.  0. CO )      PS  I      =     PI     -     P3 I 

COMPUTE     LEFT-HAN!  C     A'-'EA     LOSS 

TH'JT:     =     DMAX  1  ( -THETA,      ^SI-A^^C) 
TH-T2     =     DVAXl ( -T HTTP,      TH3T1) 
IF      (  (  THLT2-THF  Tl  )  .L  r.PTOL)      ,',  C     T2     171 
AOJT     =      ((ThET2     -    THETl)*     pmx**2    -     .J*?2 
X        -     D7AN( THFTl ) ) )     /     2.00 


*      (■)TAN(ThET2) 


ccmput::  central  area  lcss 

THFTl     =      THET2 

THFT2     =     OMAXlCTHcTl,      TH^ITP) 

THET2     =     0MIN1(THET2i      PSI+A9C) 

AOJT     =     AOUT  ♦■(T  l-i^T2-TI-ET  I  )  »(KMX<-f^MN  )  *  (PMX-RVN  ) /- 

COMPUTE     PIGMT-HANO     LU3S 

THZTl     =     THET2 

THZT2    =      OMA  Xl(  Tf   ET  1  ,      THE  TA  ) 
THET2     =      OMINl   (  Th-ET2  ,      PSI+Ar;C) 
IF      ((THETP-THETl).Lr.KTOL)      GC     TO     175 
AQJT  =  A.TUT  +  (  (  THET  2-THtT  !  )  *{  r?M  X**2)-(  D**  2)  *  (  r.TANC 
X  -     OT  AN(TfiET  1  )    )  )     /     2.00 

CONT  IN'jr 

IF      (FACT  OP)      17  6.      17-^.      175 
VrjL2     =     VnL2     -      (VCL     -     AOUT  ) 
GT     TO     2  00 

VOL2     =     VCL2     -     AOUT 
C  TITINUF 
IF     (  VOLl  .LT  .r^TOL  )     GfJ     TO     800 

CM:^     of     CHdCK     FC:^     OOUNCAnY     POSITION 


.00 


ThETZ  ) 


OUT      =    FALSE 

W  T     =      ,v  T     *     VOL  2 

I  P  T     =      I P  T     +     1 

WRl  Tf      (  -^.2  000  0       IPT 


XI ,     Y 1 ,     Zl 


mT 


C^INT  INiJE 

Cf  NT1\UF 

IF      (LiUT   )      GO     TO     q.  CO 

CONT  I  WJC 

C'.J-IT  I  NUF 

W  I  T  t":     (  ^.  .  1  0  0  0  0  )       I  f '  T 

'"•'ETUKN 
JO     f^O   !M  AT(   1  M-,    •  TO  TAL     fi'n'f^v.f       CF      IMTOGPATICN     '^C  I  fj  I  S  : 
0  0      hM-?MAT(lM      ,  15  .  IP -502  .;.:/)  ) 
0  0     -~0  ->  M  A  T  (   ^^  1  ,  '      P  N  r      •   ,  i  0  X  ,  •  X  ♦    ,  1  SX  ,  •  Y  '    ,  1  <:  X  ,  '   .""    .  1  '-  X  , 

S'iMRCUT  I^JZ     INTG-TCZ.     N3KIP.      P.      WX,      .v  T  Z  .     rsp  ^ ,      £C 


'   .13) 
•  .V  T  •  /  ) 
ALU  ) 
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Table   D-1 .    (Continued) 


C 

c 
c 
c 
c 


c 

C 


1  00 
200 


400 


I'lr^LICIT     REAL*  yC  A-H,  C-Z) 

THJ^S     SjyRQUTINE     C-ENTRATES     A     TA^JLAP     GP  10     OF     P  C  I  N  7  S 
U5-.J     FOP     ThL     NUMERICAL      INTEGPATION     PPCC=='DUrE     '.-J  I T  F 
EVFTRY     NSKIP     POINTS     OF     THE     HE  R"-^  A\-3K  I  LL  ^'A  N     VEGH. 

OI'^ENSIGN     P(l),      PX(1),      v.TZ(l) 
DAT^     '^I2/'7.gec6     C4401      03936     DOO/ 
RTRAN(XXX  ,SSS)      =     -DLGG ( X XX  )/ SSS 

SCALE     =     {  9.00<:PI  2/(  1  .  2302*2)  )**  (  3.  33333^3333^^^3^r-01  ) 
N^T     -     1  ■  ^-^      Ji  / 

RXCJPT)      =    0.00 

P  (  M  '^  T  )      -      1  .  D  0 

X     -     0.00 

OX     =     2. 50-3     *      SCALE 

SC^LE    -      1  .DO     /     (  DX     «     0rL0AT(441)) 

IJ     =     0 

GENERATE     EVERY    N  SK  I  ^     PJINT     GF     HE  RM  A.NI-S  KI  LL  ^  A  N     GRIO 

DO     200     1-1,10 
00     10  0     J  =  1  , 4  0 
X    =     X     +     OX 
IJ     =     IJ     ♦•      1 

IF      (,V00(  I  J.NSK  IP  ).NE.C  )      GC     TO     100 

NOT     =     MPT     •^     1 

■?X(^JPT)      -      X 

R(MPT)     =        DEXP  (-SCALe-*X) 

CONTI NUE 

OX    =     2.0  0     *    CX 

RX (NPT+1 )      =     X 

P(NJPT<-1)      =     COO 

^O.VM     =     — ?TRAN((1.D0     +     R  (  2  )  )  /  2.  O  0  ,      SC  Al   E  ) 

O  0     '1. 0  0     I  =  1  ,  N  P  T 

RJD     =     RTPAN((R(1)     +     R( I+l ) )/2.O0,      SCALE) 

V;T  ZC  I  )     =     RUP     -     .-^CWN 

RO'VN     =     RUP 

RETURN 

ENO 
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